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Abstract 

AFIT  research  in  support  of  the  Advanced  Logistics  Project  is  directed  at 
developing  a  Mission-Resource  Value  Assessment  Tool  for  rationally  assigning  relative 
value  to  resources  and  identifying  alternative  force  mixes  to  logistics  and  operational 
planners.  Research  of  factors  that  affect  force  mix  composition  has  been  strictly  limited 
to  how  the  operating  environment  of  USAF  combat  aircraft  influences  their  performance 
in  specified  aerospace  missions.  In  contrast,  this  research  makes  use  of  an  aircraft's 
designed  suitability  to  perform  specified  aerospace  missions  in  order  to  examine  the 
tradeoff  between  mission  suitability  and  the  amount  of  lift  needed  to  deploy  and  operate 
the  asset. 

An  Evolutionary  approach  was  applied  to  a  tri-objective  constrained  optimization 
problem  with  1 5  decision  variables  with  the  goal  of  producing  five  Pareto  optimal  sets  of 
force  mixes  corresponding  to  five  progressively  larger  sortie  capability  levels.  Analysis 
of  the  results  includes  absolute  performance  comparisons  using  different  operating 
parameter  settings,  and  time  complexity  in  relation  to  problem  scale.  Preliminary  results 
were  also  generated  from  a  version  of  the  algorithm  that  uses  a  solution  repair  function. 
These  results  help  to  assess  the  viability  of  using  a  multi-objective  fast  messy  genetic 
algorithm  to  identify  well  balanced  force  mixes. 


Know  your  limits.  Now  shatter  them. 
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IDENTIFICATION  OF  PREFERRED  OPERATIONAL  PLAN 


FORCE  MIXES  USING  A  MULTIOBJECTIVE  METHODOLOGY 
TO  OPTIMIZE  RESOURCE  SUITABILITY  AND  LIFT  COST 


I.  Introduction 


Background 

Death  by  inventory  is  the  concept  of  stockpiling  excessive  inventory  to 
compensate  for  poor  logistics  management.  This  is  a  difficult  and  expensive  business 
practice,  and  unfortunately,  is  employed  in  DoD  operations.  Take,  for  example,  the 
strategic  movement  of  personnel  and  equipment  during  the  Gulf  War.  Although  branded 
a  success  by  our  leaders  in  DoD,  the  war  served  to  highlight  logistical  problems.  “One  of 
USTRANSCOM’s  most  intractable  and  high-visibility  problems  during  Desert 
Shield/Desert  Storm  was  a  backlog  of  sustainment  cargo  at  aerial  pots  of  embarkation, 
primarily  in  the  United  States”  (Matthews  and  Holt,  1996:84).  The  amount  of  arriving 
cargo  also  overwhelmed  the  destination  points.  “Half  of  the  40,000  bulk  containers 
shipped  into  the  theater  had  to  be  opened  in  order  to  identify  their  contents,  and  most  of  it 
failed  to  contribute  in  any  way  to  our  success  on  the  battlefield”  (Muczyk,  1997:89). 
These  problems  illustrate  a  serious  gap  in  what  the  combatant  commander  or  operator 
wants  to  accomplish  and  what  the  logistician  can  make  available.  This  disconnect 
between  operations  and  logistics  is  important,  as  Paul  Judge  notes: 

The  warfighting  commander  demands  visibility  of  assets  and  requires  confidence 

in  rapid  availability.  Without  direct  knowledge  that  commodities  and  reparables 
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are  available  and  capable  of  supply  in  a  specified  time  period,  the  field 
commander  is  forced  to  stock-pile  anticipatory  requirements.  (1998:25) 

Judge  also  noted  that  recent  studies  concerning  modem  defense  logistics  supports 

a  solution  that  implements  “a  near  real-time  information  system  that  cuts  across  all 

logistic  functional  areas  and  is  not  excessively  dependent  on  manual  entry  for  raw 

information”  (1998:8). 

The  Advanced  Logistics  Project  (ALP)  is  an  effort  by  the  Defense  Advanced 
Research  Projects  Agency  (DARPA)  to  develop  a  distributed  computing  architecture  that 
links  current  and  planned  logistics  information  systems  to  the  deliberate  and  crisis  action 
planning  processes,  databases,  and  policies— an  end-to-end  system  linking  operations  and 
logistics. 
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Figure  1 .  ALP  Operational  Concept  (Carrico,  2000) 
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The  deliverable  program  is  expected  to  give  logistics  planners  a  tool  that  uses 
real-time  data  to  rapidly  develop  a  campaign  specific  logistics  plan  and  perform  dynamic 
replanning  (Carrico,  2000). 

ALP  leverages  the  revolution  in  technology  to  narrow  the  gap  between  operations 
and  logistics,  giving  planners  the  capability  to  review  multiple  deployment  plans  that  are 
based  on  current  information.  But  when  comparing  numerous  plans  and  scenarios,  how 
does  one  answer,  “What  is  best?”  In  other  words,  from  a  pool  of  available  resources,  is 
there  a  mix  of  those  resources  that  would  be  preferred  by  the  combatant  commander  over 
all  others?  AFIT  research  in  support  of  ALP  is  directed  at  developing  a  Mission- 
Resource  Value  Assessment  Tool  (M-R  VAT)  that  “rationally  assigns  relative  value  to 
material  resources”  (Swartz,  1999)  and  identifies  alternative  force  mix  compositions  to 
planners. 

A  Mission  Ready  Resource  (MRR)  is  a  combination  of  an  asset  type  and  its 
resources,  e.g.  aircraft,  pilot,  fuel,  munitions,  support  equipment  and  personnel,  etc.,  that 
is  designed  to  have  a  certain  suitability  for  a  single  task.  A  combination  of  MRR  types  is 
defined  to  be  a  MRR  set  or  force  mix.  To  demonstrate,  assume  that  a  notional  aircraft  F, 
has  two  configurations,  FA  and  FB,  which  constitutes  two  MRR  types.  Further,  if  the 
aircraft  could  be  prepared  and  flown  three  times  per  day,  then  it  would  represent  three 
MRRs  per  day.  These  three  MRR’s  could  either  be  all  FA  configurations,  or  all  FB 
configurations,  or  some  combination  of  the  two  configurations.  A  MRR  is  consumed  in 
the  performance  of  its  task,  i.e.  fuel,  munitions,  engine  cycles,  etc.  The  goal  of  a 
logistics  plan  is  to  provide  the  combatant  commander  with  the  Mission  Ready  Resource 
(MRR)  sets  that  satisfy  the  time-phased  need  for  those  resources.  In  other  words,  the 
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combatant  commander  desires  to  fulfill  a  task  with  the  best  combination  of  MRRs  that  is 


dependent  on  resource  suitability  to  a  given  task,  resource  level,  and  time,  along  with 
theater-specific  factors.  There  is  most  likely  a  number  of  acceptable  MRR  sets  for  which 
certain  combinations  of  MRRs  may  be  assessed  as  having  equivalent  suitability.  It  is 
desirable  for  the  planner  to  select  the  MRR  combination  that  best  maximizes  the  time- 
and  resource  level  dependent  suitability. 

An  MRR  set  provides  a  certain  suitability  and  capability  to  the  combatant 
commander,  but  at  a  cost:  consumption  of  lift  resources.  The  finite  lift  capability  of  the 
U.S.  military  is  a  key  constraint  on  the  amount  and  timing  of  resources  flowing  into  a 
target  area.  Therefore,  it  is  desirable  to  deliver  an  MRR  combination  that  minimizes  lift. 
There  are  then  two  competing  objectives  that  the  planner  must  deal  with  in  order  to 
present  the  best  MRR  combination  to  the  decision  maker:  maximize  asset  suitability  and 
minimize  lift  consumption. 

Problem  Statement 

Given  a  choice  among  time-phased  asset  sets,  simultaneously  minimize  lift 
resource  consumption  (cost)  and  maximize  asset  set  suitability  over  time. 

Research  Questions 

1 .  Which  methodologies  can  be  used  to  trade-off  lift  cost  and  asset  suitability  and  result 
in  an  asset  mix  that  is  preferred  over  others? 

2.  What  are  the  forms  of  decision  and  objective  spaces? 

3.  How  should  the  multiobjective  optimization  approach  be  evaluated? 
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4.  Does  the  selected  approach  result  in  an  acceptable  solution  to  the  research  problem  in 
a  reasonable  amount  of  time? 

Research  Methodology 

The  two  objectives,  maximize  suitability/minimize  lift,  are  in  conflict,  so  it  may 
be  that  no  single  optimum  solution  exists  with  respect  to  both  objectives.  It  is  desirable 
to  ascertain  whether  an  exact  solution  can  be  obtained  within  reasonable  time.  If  not,  an 
acceptable  compromise  solution  must  be  found.  Candidate  approaches  exist  and  are  well 
documented  in  literature.  Candidate  approach  selection  is  based  on  its  applicability  to  the 
problem,  its  overall  utility,  and  its  utilization  of  computational  resources. 

For  any  acceptable  asset  set,  the  combination  of  lift  costs  for  all  asset  sets  over  a 
given  period  is  constrained  by  the  maximum  throughput  of  the  transportation  pipeline 
during  that  period.  For  the  initial  model  formulation,  it  is  assumed  that  the  combined  lift 
will  not  exceed  the  maximum  available  lift.  A  constraint  can  be  added  to  subsequent 
models  in  a  way  that  limits  asset  set  selection  to  allocated  lift  consumption. 

Assumptions 

The  combinatorial  nature  of  this  research  necessitates  the  use  of  a  relatively  small 
set  of  assets  with  which  to  explore  the  algorithmic  search  for  an  acceptable  solution. 
Additionally,  actual  asset  capability  with  regard  to  specific  missions  may  be  classified. 
Therefore,  this  research  is  constructed  around  a  notional  Air  Expeditionary  Force, 
comparable  in  size  and  diversity  to  that  depicted  in  the  ALP  Pilot  Problem  (Swartz, 

1999).  No  actual  assessments  of  aircraft  suitability  or  lift  cost  will  be  used.  The  results 
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of  this  research  may  be  sensitive  to  scale  in  terms  of  time  required  to  produce  a  set  of 
acceptable  solutions. 

It  is  difficult  to  determine  what  the  actual  logistical  footprint  is  for  a  given  asset 
set.  The  logistics  support  for  a  force  is  often  undefined  until  just  hours  prior  to 
movement  (Judge,  1998:32).  It  was  noted  during  Desert  Storm  that  “the  actual  material 
shipped  grew  in  size  without  anyone’s  knowledge  and  certainly  without  any  tools  to 
predict  the  eventual  impact”  (Lynn,  1997:15).  To  be  of  use,  this  research  assumes  that 
asset  set  lift  consumption  is  both  sufficiently  accurate  and  available  for  planning. 
Research  conducted  by  Matt  Goddard  suggests  that,  for  F-16s,  the  relationship  of  asset 
quantity  to  consumption  is  linear  (2001). 

Scope/Limitations 

In  their  research  on  a  campaign  planning  decision  support  tool,  Christopher  Buzo 
and  Paul  Filcek  proposed  a  comparison  of  competing  sets  of  combat  aircraft  assets  based 
on  two  criteria  (Buzo,  2000:2,  Filcek,  2001).  The  first  criterion  is  the  intrinsic  suitability 
of  a  MRR  to  one  or  more  specific  missions.  For  instance,  a  properly  configured  F-16C  is 
capable  of  adequately  performing  many  roles  such  as  Air  Interception,  Short-Range 
Reconnaissance,  and  Air-to-Air.  In  contrast,  a  B-2’s  intrinsic  suitability  is  relatively 
confined  to  Strategic  Bombing. 

The  second  criterion  is  based  upon  situation  specific  or  extrinsic  factors  of  the 
campaign  itself.  Such  factors  modify  the  task  suitability  of  MRRs  in  certain  scenarios. 
For  example,  a  planner  with  knowledge  of  an  enemy  with  a  high  anti-air  capability  would 
alter  a  force  mix  in  favor  of  aircraft  with  a  high  absolute  suitability  in  air  defense 
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suppression.  Political  issues  are  also  extrinsic  factors  of  a  campaign,  forcing  planners  to 
consider  over-flight  rules,  coalition  participation,  and  beddown  constraints,  to  name  a 
few. 

Buzo’s  and  Filcek’s  research  of  factors  that  would  affect  force  mixes  was  strictly 
limited  to  extrinsic  factors  of  USAF  combat  aircraft  to  perform  specific  aerospace 
missions  (2000,  2001).  In  contrast,  this  research  makes  use  of  the  intrinsic  suitability  of 
MRRs  to  examine  the  tradeoff  of  between  MRR  suitability  against  MRR  lift  cost — 
extrinsic  factors  are  not  considered. 

Discussion  of  operational  plans  and  actual  task  suitability  may  be  classified. 
Therefore,  this  research  and  its  conclusion  is  restricted  to  the  unclassified  realm. 

Summary 

This  chapter  illustrated  the  disconnect  between  operations  planners  and  logistics 
planners  and  ALP’s  intention  to  close  that  gap.  AFIT  research  is  centered  on  providing  a 
decision  support  tool  that  would  help  campaign  planners  select  the  best  force  mix  from  a 
pool  of  combinations.  The  problem  of  force  mix  selection  can  be  handled  as  a 
multiobjective  optimization  of  two  competing  objectives:  minimize  lift  resource 
consumption  (cost)  and  maximize  MRR  task  suitability.  This  research  proposes  to 
develop  a  methodology  for  presenting  the  decision  maker  with  a  set  of  acceptable  force 
mixes  after  which  an  extrinsic  assessment  is  then  made  and  the  user-defined  best  force 
mix  is  selected. 

Chapter  II  provides  a  background  on  multiobjective  optimization,  discusses  global 
versus  local  optimization,  and  reviews  classic  and  modem  multiobjective  techniques. 
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Chapter  III  describes  the  methodology  used  to  construct  the  multiobjective  problem  and 
evaluate  the  solution  approach.  Chapter  IV  details  the  results  using  the  selected 
multiobjective  optimization  methodology.  Chapter  V  provides  conclusions  on  research 
contributions  and  makes  recommendations  for  further  research. 


8 


II.  Literature  Review 


Introduction 

Mathematical  optimization  is  the  formal  title  given  to  the  branch  of  computational 
science  that  seeks  to  answer  the  question  ‘What  is  best?’  for  problems  in  which  the 
quality  of  any  answer  can  be  expressed  as  a  numerical  value.  Such  problems  arise  in  all 
areas  of  mathematics,  the  physical,  chemical,  and  biological  sciences,  engineering, 
architecture,  economics,  and  management,  and  the  range  of  techniques  available  to  solve 
them  is  nearly  as  wide.  In  practical  problems,  we  often  want  to  optimize  more  than  one 
measure  of  performance  at  once. 

This  research  is  concerned  with  a  particular  problem  class:  multiobjective 
optimization  problems  (MOPs).  The  purpose  of  this  chapter  is  to  provide  a  broad 
overview  of  the  field  in  order  answer  the  first  research  question:  “which  methodologies 
can  be  used  to  trade-off  lift  cost  and  asset  suitability  and  result  in  an  asset  mix  that  is 
preferred  over  others?”  Following  an  overview  on  the  MOP  class,  modem  approaches  to 
handling  MOPs  are  presented. 

MOP  Overview 

The  goal  of  an  optimization  problem  can  be  formulated  as  follows:  find  the 
combination  of  parameters  (independent  variables)  that  optimize  (maximize  or  minimize) 
a  given  quantity,  possibly  subject  to  some  restrictions  on  the  allowed  parameter  ranges. 
The  objective  function  is  the  quantity  to  be  optimized.  The  parameters  that  can  be 
changed  are  the  decision  variables  that  represent  discrete  solutions  of  a  combinatorial 
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problem  (Reeves,  1995:2).  The  restrictions  on  parameter  values  are  termed  constraints. 
For  convenience  of  mathematical  treatment,  all  problems  in  this  thesis  are  assumed  to 
have  minimization  objectives  unless  stated  otherwise. 

To  solve  a  decision-making  problem  analytically,  it  is  helpful  to  state  the  problem 
in  numerical  terms.  Given  that  an  objective,  O,  has  a  corresponding  ^-dimensional  set  of 
alternatives,  X,  the  criterion  for  the  objective  is  defined  as  an  objective  function: 

R1  (2.1) 

where /is  a  mapping  that  may  be  linear  or  non-linear  (Van  Veldhuizen,  1999:2-2).  The 
general  form  of  an  MOP  with  p  objectives  is 
minimize: 

f(x)  =  (f(x),f2(x),...,fp(x))  over xeX  (2.2) 

subject  to: 


c,00  =  o 

i  = 

(2.3) 

c,.(x)>0 

i  =  m  +\,...,m 

(2.4) 

where  x  is  the  column  vector  of  the  n  independent  variables,  and  c,(. x)  is  the  set  of 
constraint  functions  that,  depending  on  the  situation,  may  or  may  not  be  included  in  the 
problem  (Sawaragi,  et  al.,  1985:2). 

In  an  MOP,  it  is  difficult  to  obtain  a  unique  optimal  solution.  This  is  because  the 
problem’s  objectives  are  usually  in  conflict  with  one  another:  one  cannot  improve  the 
performance  of  a  particular  objective  without  causing  a  corresponding  deterioration  in 
performance  in  one  or  all  of  the  others.  Examples  of  conflicting  objectives  may  be 
maximizing  speed  and  safety  in  a  vehicle,  or  minimizing  acquisition  cost  and  schedule  of 
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a  new  aircraft  while  maximizing  its  performance.  Classical  methods  of  dealing  with  this 
problem  produce  a  solution  by  combining  objectives  in  some  way  that  is  usually  a 
subjective  expression  of  an  a  priori  not  well  understood  trade-off  surface  (Fonseca  and 
Fleming,  1993:1,  Horn  et  al.  1993:1,  Deb,  1999b:4,  Ehrgott  and  Gandibleux,  2000:12, 
Taber  et  al.,  1999:1).  By  treating  the  problem  as  multiple,  competing  objectives,  the 
result  is  a  set  of  solutions  in  decision  variable  space  whose  components  represent  a  trade¬ 
off  in  the  objective  function  space.  A  decision  maker  implicitly  chooses  an  acceptable 
solution  (or  solutions)  by  selecting  one  or  more  of  these  alternatives  (Van  Veldhuizen, 
1999:2-2).  Ideally,  the  alternatives  should  be  selected  from  a  set  of  equally  preferred 
solutions  called  the  Pareto  optimal  set,  or  simply  P.  Pareto  optimal  solutions  are  also 
termed  efficient  or  admissible  solutions  (Yu,  1985:22).  The  mapping  of  P  to  the 
objective  space  forms  the  Pareto  front,  PF.  The  Pareto  front  is  also  known  as  the 
nondominated  front.  The  reader  is  referred  to  Appendix  A  for  additional  discussion  of 
Pareto  concepts. 

Global  Optimization.  The  desired  solution  to  an  optimization  problem  is  the  true 
or  global  optimum.  The  strict  definition  of  the  global  optimum  (minimum)  x'  of f{x)  is 

/ (*')  <  f(y )  for  all  y  e  V (x),  y  *  x',  (2.5) 

where  V(x)  is  the  set  of  feasible  values  of  the  decision  variables  x  (Allen,  et  al.,  1 996).  A 
complication  that  arises  in  nonlinear  optimization  is  that  a  local  optimum  need  not  be  a 
global  optimum.  For  example,  consider  the  function  of  a  single  variable  plotted  in  Figure 
2.  Over  the  interval  0  <  x  <  20,  this  function  has  three  local  maxima — x  =  3.8,  x  =  10.2, 
and  x  =  16.7 — but  only  one  of  these — x  =  3.8 — is  the  global  maximum. 
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Figure  2.  Global  and  Local  Maxima 

Finding  the  global  optimum  of  a  general  mixed  integer  MOP  is  NP-Complete  and 
MOP  solutions  that  satisfy  all  constraints  and  globally  optimize  all  objective  functions 
may  not  even  exist  (Van  Veldhuizen,  1999:2-2). 

Techniques  for  Solving  MOPs.  Techniques  for  solving  MOPs  have  existed  for 
years  and  usually  rely  upon  either  enumerative  or  approximation  approaches  (Ehrgott  and 
Gandibleux,  2000:12,  20).  For  MOPs,  it  is  from  the  Pareto  optimal  set  that  an  informed 
decision  maker  chooses  a  compromise  solution.  Ideally,  the  presented  set  is  the  true 
Pareto  optimal  set  P,rue  that  corresponds  with  the  true  Pareto  front  PFtrUe-  However, 
complete  enumeration  of  solutions  for  even  a  reasonable  MOP  is  impractical  from  both 
computational  and  decision  making  standpoints. 

A  MOP  is  required  to  pare  down  the  solution  set  to  one  in  which  the  decision 
maker  can  feasibly  use  to  select  a  solution  that  represents  the  best  tradeoff  between 
objectives.  In  a  real-world  problem  with  real-valued  solutions,  the  presented  Pareto 
optimal  set  Pknown  is  a  discretized  approximation  or  subset  of  a  continuous  P,n,e- 
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According  to  Reeves,  this  kind  of  solution  tends  to  favor  an  approximate  or  heuristic 
approach  to  finding  it  (1995:1 1).  While  an  exact  model  to  a  real-world  problem  is 
beyond  our  reach,  “it  may  be  possible  to  model  the  real-world  problem  rather  more 
accurately  than  is  possible  than  if  an  exact  algorithm  is  used”  (Reeves,  1995:1 1).  A 
heuristic  allows  us  to  solve  optimization  problems  “in  ways  that  are  less  than  perfect  yet 
of  considerable  practical  value”  (Harel,  1987:344). 

Deterministic  heuristics,  whose  members  include,  but  are  not  limited  to,  greedy 
algorithms,  descent  algorithms,  and  deterministic  linear  /  non-linear  programming 
methods,  use  problem  domain  knowledge  to  shrink  the  solution  space  (Van  Veldhuizen, 
1999:2-10).  For  any  heuristic,  a  reduction  in  computational  cost  comes  without  being 
able  to  guarantee  either  feasibility  or  optimality.  Further,  deterministic  algorithms,  when 
applied  to  MOPs,  suffer  from  their  poor  handling  of  irregularities  in  the  search  space — 
high-dimensional,  discontinuous,  multimodal,  and  /  or  AP-Complete — and  can  be 
expected  to  produce  only  local  solutions  (Van  Veldhuizen,  1999:1-3,  2-12,  Reeves, 
1995:6).  MOPs  are  usually  better  handled  by  flexible  and  more  robust  heuristic 
approaches  (Reeves:  1995:1 1). 

When  systematic  search  methods  fail,  stochastic  techniques  are  used.  The 
algorithms  discussed  in  the  following  sections  all  employ  some  form  of  random  search. 

Modern  Methods  for  Handling  MOPs 

Heuristic  approaches  are  typically  designed  for  a  specific  problem  and  are  not 
suited  for  a  wide  range  of  applications  (Ehrgott  and  Gandibleux,  2000:15).  In  contrast,  a 
metaheuristic  employs  a  master  strategy  to  take  advantage  of  the  search  space  and  guide 
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the  search  (Ehrgott  and  Gandibleux,  2000: 15).  Metaheuristics  are  much  more  general  in 
their  application.  Recent  advances  in  computational  power  have  pushed  metaheuristics  to 
the  forefront.  Simulated  Annealing  (SA),  Tabu  Search  (TS),  and  Evolutionary 
Algorithms  (EAs)  are  examples  of  metaheuristics  that  have  been  researched  in  great 
depth  and  have  produced  superior  results,  in  terms  of  both  solution  quality  and 
computational  efficiency,  in  a  wide  variety  of  applications  (Ehrgott  and  Gandibleux, 
2000:15).  All  methods  rely  heavily  on  computing  skills  for  practical  implementation. 

There  are  two  main  approaches  used  by  metaheuristics:  1)  local  search  in 
objective  space  and  2)  population  based. 

Local  Search  in  Objective  Space.  Based  on  the  principle  of  search  directions,  this 
approach  starts  from  some  initial  solution  and  proceeds  in  a  given  search  direction  to 
focus  on  a  portion  of  the  nondominated  front.  The  search  proceeds  iteratively  in  other 
search  directions  in  order  to  approximate  the  entire  Pareto  front.  “At  any  time  the  search 
mechanism  uses  only  one  solution  and  an  iteration  tries  to  attract  the  solution  generated” 
towards  the  Pareto  front  along  the  given  direction  (Ehrgott  and  Gandibleux,  2000:16). 

Hill  Climbing,  Simulated  Annealing  (SA)  and  Tabu  Search  are  examples  of  the 
first  approach.  Hill  climbing  begins  with  a  single  random  solution  that  is  perturbed  to 
change  its  evaluation.  After  several  such  perturbations,  the  best  evaluation  is  chosen  as 
the  next  starting  point.  Continuing  on  in  this  way  eventually  results  in  reaching  an 
optimum.  However,  it  is  not  known  whether  it  has  reached  local  or  global  optimum.  The 
search  space  can  be  explored  by  starting  repeatedly  with  a  new  random  solution  in  hopes 
of  finding  a  better  solution  (Caryl,  no  date). 
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SA  employs  an  analogy  between  the  way  in  which  a  metal  cools  and  freezes  into 
a  minimum  energy  crystalline  structure,  and  the  search  for  a  minimum  in  a  general 
system.  At  high  temperatures,  the  molecules  of  a  liquid  move  freely.  As  the  liquid  cools, 
that  mobility  becomes  restricted  and  the  molecules  achieve  crystalline  form  and  the 
system’s  minimum  energy  state,  i.e.  global  optimum. 

Like  hill  climbing,  SA  randomly  perturbs  the  objective  function  in  a  way  that 
will,  for  a  given  change,  causes  a  decrease  and  for  another  change  causes  an  increase. 

But  instead  of  selecting  the  best  evaluation  and  continuing  on  (analogous  to  rapidly 
cooling  the  liquid),  SA  introduces  an  element  of  randomness:  a  change  that  does  not 
improve  upon  the  current  optimum  is  executed  with  a  probability  p  <  1 .  This  is  typically 
based  on  the  Boltzmann  probability  distribution: 

E 

P(E )  ~  ekT 

where  E  is  the  energy  of  the  state,  k  is  Boltzmann’s  constant,  and  T  is  the  temperature. 
This  equation  means  that  the  probability  of  finding  a  particle  with  energy  E  is 
proportional  to  the  exponential  of  -E  divided  by  the  product  of  k  and  T  (Rappe,  no  date). 
So  at  a  given  temperature,  a  system  can  be  in  a  range  of  possible  energy  states.  A  higher 
temperature  increases  the  likelihood  of  a  high  energy  state.  Simulated  annealing  makes 
use  of  the  fact  that  at  low  temperatures,  there  is  still  of  chance  of  being  in  a  high  energy 
state,  thus  allowing  a  jump  out  of  a  local  minima.  This  random  perturbation  gives  an  SA 
its  ability  to  avoid  being  trapped  at  a  local  optimum  (Reeves,  1995:26). 

An  SA  implementation  tends  to  be  problem  specific.  The  choice  of  the 
temperature  or  annealing  schedule  depends  on  the  expected  range  of  function  values  and 
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the  shape  of  the  function  surface.  Experimentation  is  required  to  obtain  a  method  that 
works  well  for  a  particular  problem. 

TS  is  can  also  be  employed  as  a  stochastic  method,  but  rather  than  using  random 
moves,  TS  employs  a  directed  search  along  with  a  memory  to  imitate  intelligent 
processes  (Reeves,  1995:13).  TS  is  a  form  of  neighborhood  search.  Beginning  with  a 
solution  within  a  defined  neighborhood,  the  algorithm  proceeds  iteratively  to  visit  a  series 
of  locally  optimal  solutions.  At  each  iteration,  a  best  neighbor  is  chosen  to  replace  the 
current  solution.  To  allow  the  search  to  move  beyond  local  optima,  a  list  of  moves  that 
are  not  allowed  or  tabu  is  used.  This  list  prevents  recently  visited  solutions  from  being 
considered  for  a  given  number  of  iterations  of  the  algorithm  (Reeves,  1995:83,  86-88). 

Population  Based  Approach.  The  population  based  approach  takes  advantage  of 
information  carried  by  a  population  of  solutions.  Heuristics  using  this  method 
predominantly  fall  under  the  category  of  Evolutionary  Algorithms  (EAs),  a  class  that  uses 
the  evolutionary  concepts  of  survival  of  the  fittest  and  generational  improvement  as  its 
inspiration. 

Despite  the  probabilistic  nature  of  EA  operators,  EAs  are  not  completely  random 
searches  and  are  directed  by  the  information  carried  by  the  population.  Contrary  to  the 
local  search  methods,  where  only  one  individual  is  attracted  toward  the  Pareto  front,  here 
the  entire  population  contributes  to  the  evolutionary  process  toward  the  Pareto  front  by 
searching  for  many  nondominated  solutions  at  once.  It  is  this  characteristic  that  makes  the 
population  based  approach  very  attractive  for  solving  MOPs,  but  it  comes  at  a  relatively 
higher  computational  cost  since  many  fitness  evaluations  are  required  (Ehrgott  and 
Gandibleux,  2000:16,  Fonseca  and  Fleming,  1993:1). 
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Recent  research  by  Van  Veldhuizen  catalogued  206  multiobjective  EAs  (MOEAs) 
(1999:A-1).  A  widespread  implementation  of  the  population-based  approach  is  a  type  of 
EA  called  the  Genetic  Algorithm  (GA).  Its  broad  applicability,  ease  of  use,  and  success 
in  handling  MOPs  makes  it  “no  surprise  that  a  number  of  different  multi-objective  GA 
implementations  exist  in  the  literature. . (Deb,  1 998:4).  A  GA,  in  a  single  run,  can 
provide  “a  number  of  Pareto-optimal  solutions”  (Deb,  1998:26)  and  has  the  “ability  to 
find  global  optima  while  being  able  to  cope  with  discontinuous  and  noisy  functions” 
(Fonseca  and  Fleming,  1993:7).  A  detailed  discussion  of  GA  theory,  operation,  and  types 
of  MOEAs  can  be  found  in  the  next  section. 

The  limit  of  any  optimization  method  is  succinctly  expressed  by  the  No  Free 
Lunch  theorems  (Wolpert  and  Macready,  1996),  which  tell  us  that  “used  blindly,  there  is 
an  equal  chance  that  any  optimization  technique  will  perform  the  same”  (Practical  Guide, 
no  date).  The  choice  of  which  optimization  method  to  use  should  be  based  on  what  is 
known  about  the  system  being  optimized.  Even  though  a  GA  has  no  guarantee  of 
performing  better  than  another  method  in  a  given  application,  in  most  cases  a  GA’s 
parameters  and  configuration  can  be  tailored  to  achieve  adequate  search  performance. 

Genetic  Algorithms 

The  Genetic  Algorithm  (GA)  is  based  on  the  biological  processes  of  evolution  as 
described  by  Charles  Darwin  (1936).  An  organism  is  made  up  of  genetic  material 
embedded  with  environmental  knowledge.  Natural  selection  sees  to  it  that  those 
individuals  that  are  better  suited  to  their  environment  survive  to  pass-on  their  good 
genetic  material.  The  less  fortunate  die  and  take  their  bad  genes  with  them. 
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Reproduction  happens  in  an  environment  where  the  selection  of  who  gets  to  mate  is 
largely  a  function  of  individual  fitness.  Reproducing  pairs  or  parents,  produce  offspring 
with  chromosomes  containing  information  from  each  parent.  Evolution  uses  mutation  to 
stimulate  diversity  in  the  population.  Mutated  individuals  do  not  always  survive,  but 
occasionally  there  are  those  that  are  better  suited  to  their  environment  and  more 
competitive  than  the  others.  Their  environmental  advantage  is  passed  on  to  their 
offspring  and  ultimately  to  future  generations  (Caryl,  no  date). 

Basic  Operation.  There  are  five  characteristic  components  in  every  GA  (Caryl): 

1 .  A  way  to  create  an  initial  population  of  potential  solutions 

2.  A  genetic  representation  for  solutions  to  the  problem 

3.  An  evaluation  function  that  plays  the  role  of  the  environment, 
rating  solutions  in  terms  of  their  fitness 

4.  Genetic  operators  that  alter  the  composition  of  children  during 
reproduction 

5.  Values  for  various  parameters  that  the  genetic  algorithm  uses 
(population  size,  probability  of  applying  genetic  operators) 


The  pseudo  code  for  the  basic  algorithm  is  presented  in  Figure  3: 

Start  GA 

//  start  with  an  initial  time 
t  =  0; 

//  initialize  random  population  of  individuals 
initialize  P(t ); 

//evaluate  fitness  of  all  initial  individuals  of population 
evaluate  P(t ); 

//  test  for  termination  criterion 
while  not  done  do 

//  increment  the  time  counter 
t  —  t  +  1 ; 

//select  a  sub-population  for  offspring  production 
P'{t)  =  parents  from  P(t); 


18 


//recombine  the  "genes"  of  selected  parents 
recombine  P'(t ); 

//perturb  the  mated  population  stochastically 
mutate  P'{t)\ 

//  evaluate  its  new  fitness 
evaluate  P'(t ); 

//select  the  survivors  based  on  fitness 
P{t)  =  survive  P(t),  P'(t ); 

end  do 
end  GA. 

Figure  3.  Pseudo  Code  for  a  Simple  Genetic  Algorithm  (Osyczka  and  Kundu,  1995:95) 


A  single  iteration  of  the  while  loop  constitutes  a  generation. 

Initialization.  A  GA  maintains  a  population  of  solutions,  all  of  which  are 
potential  parents.  The  first  generation’s  population  is  initialized,  usually  with  randomly 
generated  individuals.  Another  technique  is  to  initialize  the  population  with  high-quality 
solutions.  This  approach  has  shown  to  increase  the  speed  of  convergence,  but  with  an 
increased  possibility  of  premature  convergence  to  only  a  portion  of  the  Pareto  front 
(Reeves,  1995:164). 

Representation.  A  GA  can  be  thought  of  as  a  “DNA  simulator”  where 
nonbiological  structures  could  be  modeled  in  terms  of  bit  strings  that  could  be  changed 
through  transformations  analogous  to  evolution.  In  the  traditional  GA,  a  genotypic 
representation  scheme  encodes  n  decision  variables  into  n  sequences  of  binary  bits  that 
together  form  a  bit  string  or  chromosome  that  represents  an  individual  solution  in 
decision  space.  The  values  of  the  variables  are  termed  alleles  and  the  location  of  a  bit 
within  the  string  is  its  locus  (Reeves,  1995:138). 

Genotype  is  the  genetic  phrase  for  the  encoded  variables  while  phenotype  is  used 
for  decoded  variable  representation.  This  encoding  scheme  leads  naturally  to 
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representation  of  integer  decision  variables.  For  example,  equation  (2.6)  is  an  objective 
function  in  two  variables: 

max:  /  (x,  y)  =  x2  +  y2  (2.6) 

Each  variable  can  be  represented  by  bit  string  of  length  m.  If  the  phenotypic 
representations  for  x  and  y  are  3  and  6  respectively  (i.e.  x  =  3  and  y  -  6  and  m  =  4,  then 
(2.7)  is  the  8-bit  string  genotypic  representation  of  x  and  y: 

0  0  11  0  110  (2.7) 

As  this  also  represents  an  entire  solution  to  equation  (2.6),  the  bit  string  is  itself  a 
chromosome. 

Continuous  variables  are  approximated  through  a  scaling  function.  The  accuracy 
with  which  an  optimum  solution  can  be  resolved  depends  on  the  length  of  the  bit  string. 
For  instance,  a  variable  represented  by  a  22  bit  chromosome  can  range  between  0  and 
222  - 1 .  When  rescaled  into  real  numbers  with  a  range  of  3,  this  representation  gives 
quantization  errors  of  7.2  x  10’7.  The  scaling  operation  also  serves  to  designate  the  upper 
and  lower  bounds  of  the  decision  variable  (Osyczka  and  Kundu,  1995:  95). 

Other  encoding  methods  that  have  been  used  are  real  number  representation  and 
an  integer  representation  variation.  In  real  number  representation,  the  each  decision 
variable  is  simply  represented  by  floating  point  number,  resulting  in  a  “one  gene  /  one 
variable  relationship”  (Practical  Guide,  no  date).  Using  an  integer  instead  of  a  floating¬ 
point  number  results  in  an  integer  representation.  The  disadvantage  in  using  a  binary 
representation  is  a  result  of  the  extra  steps  needed  to  decode  the  binary  string  to  a 
floating-point  number  and  back  for  each  fitness  evaluation  (Practical  Guide,  no  date). 
However,  contrary  to  real  number  representation,  the  theoretic  aspects  of  binary 


20 


representation  have  been  thoroughly  explored,  placing  performance  evaluation  of  the  GA 
on  much  more  solid  ground  (Practical  Guide,  no  date). 

Fitness  Evaluation.  Since  the  goal  of  optimization  is  to  either  maximize  or 
minimize  the  objective  function  value,  a  measure  of  solution  fitness  is  a  function  of  the 
objective  function  value.  If  the  objective  function  in  equation  (2.6)  is  used  as  a  direct 
measure  of  solution  fitness,  then  the  chromosome  in  equation  (2.7)  can  be  decomposed 
into  x  and  y  values  and  substituted  into  the  objective  function: 

/(3,6)  =  32  +  62  =45 

If  maximizing,  the  solution  (3,  6)  would  yield  a  better  fitness  than  solutions  with  lesser 
values. 

Once  a  population  has  been  produced,  it  can  be  evaluated  using  an  objective 
function  that  characterizes  every  individual’s  performance  in  the  problem  domain.  The 
number  of  fitness  evaluations  increases  with  the  number  of  objective  functions  and 
population  size.  Fitness  evaluation  is  the  primary  source  of  GA  computational  cost. 

Genetic  Operators.  Traditional  genetic  operators  are  selection,  recombination, 
mutation,  and  evaluation.  Selection  for  reproduction  is  dependent  on  the  evaluated 
fitness  of  each  solution,  the  idea  being  that  better  solutions  have  better  representation  in 
the  population  in  order  to  improve  the  population  with  respect  to  preceding  generations. 
Selected  solutions  are  recombined  to  form  new  solutions  that  are  then  evaluated  for 
inclusion  in  the  next  generation’s  population.  Prior  to  evaluation,  a  solution  may  be 
changed  by  probabilistically  applying  a  mutation  operation  to  one  or  more  of  its  genes. 

Selection.  Solution  fitness  is  used  to  bias  the  selection  process  toward 
highly  fit  individuals  while  still  allowing  less  fit  individuals  to  reproduce.  This  has  the 
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effect  of  keeping  a  measure  of  diversity  in  the  population  thereby  making  the  search  more 
global.  Highly  fit  individuals  are  given  a  higher  probability  of  being  selected  for 
reproduction  than  individuals  with  a  lower  fitness  value.  The  average  performance  of 
individuals  can  be  expected  to  increase  since  those  individuals  with  better  fitnesses  are 
more  likely  to  be  selected  for  reproduction  and  the  lower  fitness  individuals  are 
eventually  culled  from  the  population.  Individuals  may  be  selected  more  than  once  at 
each  iteration  of  the  GA. 

There  are  a  variety  of  selection  schemes  employed  in  GAs  (Osyczka  and  Kundu, 
1995:95,  Deb,  1999a:7).  Common  methods  include  proportionate  selection,  ranking 
selection,  and  tournament  selection.  In  proportionate  selection,  the  probability  that  an 
individual  is  chosen  for  selection  is  the  individual’s  fitness  divided  by  the  sum  of  the 
current  population’s  fitnesses: 

WO-  S*""®  (2.8) 

Y^fitness(j) 

7=1 

where  i  is  the  individual  in  question  and  n  is  the  population  size.  The  new  population  of 
potential  parents  is  then  selected  by  making  n  random  draws  from  a  uniform  distribution 
(Dorigo  and  Maniezzo,  1993:7). 

In  ranking  selection,  the  population  is  sorted  from  best  to  worst  in  terms  of 
fitness.  The  number  of  copies  that  an  individual  should  receive  is  given  by  an  assignment 
function,  and  is  proportional  to  the  rank  of  an  individual  rather  than  its  absolute  fitness 
value. 

In  tournament  selection,  successive  groups  of  k  of  individuals  are  chosen  from  the 
population  and  compared.  The  best  individual  from  each  group  is  selected  as  a  parent  for 
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the  next  generation.  This  process  is  repeated  until  the  mating  pool  is  filled.  When  k  =  2, 
each  solution  will  compete  twice.  The  best  solutions  will  have  two  copies  in  the  new 
population,  the  worst  are  eliminated,  and  those  in  between  have  one.  Deb  reports  that 
“tournament  selection  has  better  convergence  and  computational  time  complexity 
properties  compared  to  any  other  reproduction  operator  that  exist  in  the  literature,  when 
used  in  isolation”  (1999a:7). 

Recombination.  The  recombination  operation  is  typically  referred  to  as  a 
crossover.  After  selection,  each  individual  has  a  probability  pc,  called  the  crossover  rate, 
of  being  chosen  for  crossover.  This  probability  is  usually  set  high,  between  0.5  and  0.9. 
Randomly  chosen  pairs  of  individuals  are  combined  to  produce  two  offspring.  Crossover 
can  be  applied  in  different  ways.  Two  of  the  most  often  used  crossover  operators  are 
single  point  and  multi-point  crossover. 

In  single  point  crossover,  a  common  point  between  two  genes  in  both  parents  is 
selected  at  random.  The  number  of  points  is  simply  the  length  of  the  chromosome,  /, 
minus  one.  The  offspring  are  created  by  concatenating  the  pre-selection  point  genes  from 
one  parent  with  the  post-selection  genes  from  the  other  parent.  For  example,  if  the 
randomly  selected  crossover  point  is  between  the  third  and  fourth  genes  of  the 
chromosome,  as  in  Figure  4,  parents  PI  and  PI  produce  the  offspring  0\  and  02 
(Reeves,  1995:  154). 


PI 

10  10010 
t 

0\ 

1011001 

P2 

0111001 

t 

02 

0  110010 

Figure  4.  Example  of  Single-point  Crossover 
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In  multi-point  crossover,  up  to  n  -  1  points  may  be  selected  a  crossover  points,  n 
being  the  population  size.  The  parents  can  then  swap  every  other  segment.  In  uniform 
crossover,  each  gene  position  is  considered  for  crossover.  However,  this  method  has  a 
high  probability  of  producing  offspring  that  are  considerably  different  from  their  parents 
and  so  pc  is  usually  set  low,  e.g.  0.01  (Practical  Guide,  no  date). 

In  practice,  “the  simple  crossover  operator  has  proved  extremely  effective” 
(Reeves,  1995: 170).  It  does  no  backtracking  or  table  lookups,  making  it  a  simple  and 
efficient  method  to  implement. 

Mutation.  Crossover  is  responsible  for  the  search  aspect  of  the  GA  and  is 
considered  the  primary  operator.  Mutation  is  responsible  for  keeping  a  measure  of 
diversity  in  the  search  and  is  considered  as  a  background  operator.  The  offspring  from 
reproduction  are  further  perturbed  by  mutation.  Each  bit  in  a  chromosome  is  changed 
with  given  probability  pm,  called  the  mutation  rate.  In  a  binary  representation  scheme, 
this  means  flipping  the  bit. 

The  mutation  operator  is  both  simple  and  powerful  by  guaranteeing  “that  every 
point  in  the  search  space  can  be  reached”  (Dorigo  and  Maniezzo,  1993:7).  It  works  by 
biasing  the  creation  of  new  solution  in  the  neighborhood  of  the  original  solution  (Deb, 
1999a:  10).  Overuse  of  the  mutation  operator  would  destroy  this  relationship  and  so  it  is 
used  sparingly,  with  pm  usually  set  to  between  0.1  and  0.001.  Both  the  Messy  GA  (mGA) 
and  Multiobjective  Messy  GA  (MOMGA)  set  pm  to  0  and  have  obtained  highly 
competitive  results  in  comparison  with  other  GA  implementations  (Van  Veldhuizen, 
1999). 
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At  this  point,  the  next  generation’s  population  must  be  filled.  The  process 
described  by  Osyczka  and  Kundu  in  Figure  3  is  only  one  example  of  how  this  may  be 
accomplished.  Following  along,  the  new  individuals  that  result  from  crossover  and 
mutation  are  evaluated  for  fitness  and  compared  with  the  parent  population.  The  next 
generation  is  made  up  of  the  individuals  from  both  generations  with  the  highest  fitnesses. 
Deb  implements  a  variation  on  filling  the  next  generation  by  using  the  n(pc)  offspring  and 
n(  1  - pc)  parents  (Deb,  1999a:9).  This  method  is  expected  to  produce  better  solutions 
since  the  higher  fitness  parents  that  are  the  result  of  the  selection  process  are  used  to  fill 
the  vacancies  left  by  the  crossover  operator  in  the  new  population.  The  n(  1  -pc)  parents 
can  be  copied  either  deterministically  or  at  random. 

These  processes  of  selection,  recombination,  mutation,  and  evaluation  are  then 
repeated  until  some  termination  criteria  is  satisfied,  e.g.  upon  reaching  a  maximum 
number  of  generations,  a  specified  fitness,  a  specified  number  of  solutions  in  the 
nondominated  solution  set,  or  after  an  elapsed  period  of  time.  Another  stopping  rule  used 
is  based  on  the  number  of  fitness  function  evaluations  performed.  The  number  of 
function  evaluations  required  to  find  the  optimal  solution  set,  within  a  given  tolerance  of 
course,  is  an  important  measure  of  algorithm  efficiency  (Van  Veldhuizen  and  Lamont, 
2000(a):  141).  This  is  discussed  in  the  section  on  multi  objective  GAs. 

Parameter  Settings.  The  parameters  most  often  cited  as  having  a  significant 
affect  on  the  performance  of  an  EA  are  population  size  («),  crossover  rate  (pc),  and 
mutation  rate  ( pm )  (Gray,  1 997,  Caryl,  no  date.  Practical  Guide,  no  date).  Despite  the 
many  papers  on  the  theoretical  and  applied  use  of  EAs,  there  are  very  few  quantitative 
methods  for  determining  the  proper  values  to  use  in  a  given  optimization  problem 
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(Practical  Guide,  no  date).  The  parameter  values  that  produce  the  most  efficient  and 
effective  results  depend  upon  the  given  problem  and  how  the  EA  is  applied,  i.e.  search 
space  topology,  representation  scheme,  selection  and  recombination  methods.  The 
parameters  may  even  vary  with  each  generation  or  between  decision  variables 
(Mathematical  Optimization,  no  date). 

Researchers  have  used  parametric  studies  to  determine  the  best  settings  for  a 
particular  problem  (Deb  and  Agrawal,  1999:3,  Practical  Guide,  no  date).  Ad  hoc 
parameter  settings  are  based  on  what  is  generally  known  about  how  their  interaction 
affects  a  GA’s  performance,  i.e.,  algorithmic  efficiency  and  the  exploration  and 
exploitation  of  the  search  space.  For  GAs,  the  most  time-consuming  task  is  fitness 
evaluation  (Van  Veldhuizen  and  Lamont,  2000(a):  142).  GA  complexity  and  efficiency 
are  generally  stated  in  terms  of  the  number  of  fitness  evaluations  performed.  Search 
space  exploration  refers  to  how  well  population  diversity  is  maintained  in  the 
nondominated  front.  Search  space  exploitation  refers  to  how  well  the  search  is  guided 
towards  the  true  Pareto  front  (Deb,  1998:4).  The  parameter  settings  used  depend  on  what 
aspect  of  GA  performance  the  researcher  is  focused  on  (Deb  and  Agrawal,  1999:2-3). 
Swinging  a  parameter’s  settings  through  a  predetermined  minimum  and  maximum  can 
give  some  picture  of  Pareto  front.  In  consideration  of  the  interaction  between  parameter 
settings,  the  number  of  runs  needed  to  do  this  is  itself  a  multiobjective  problem  and  is 
NP-complete  (Zydallis,  2001). 

It  can  be  seen  intuitively  that  the  setting  of  the  population  size  is  a  trade-off 
between  solution  diversity  and  algorithmic  efficiency.  As  n  increases,  the  diversity  of 
among  individuals  is  expected  to  increase,  thereby  decreasing  the  chance  of  premature 
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convergence  to  suboptimal  solutions  or  only  a  portion  of  the  Pareto  front.  This  does  not 
imply  that  increasing  the  population  size  automatically  improves  convergence  to  the 
Pareto  front  (Zitzler  et  al..  1999:18).  Increasing  the  population  size  does  have  a 
computation  cost.  For  a  k  objective  optimization  problem,  at  least  kn  fitness  evaluations 
are  required.  Suggested  values  for  n  are  between  25  and  100  (Practical  Guide,  no  date). 

Deb  and  Agrawal  have  shown  that  for  simple  functions,  GAs  using  both  crossover 
and  mutation  perform  better  than  either  of  them  alone,  and  suggest  the  use  of  a  large 
crossover  probability  with  a  small  mutation  probability  (1999:20).  With  more  difficult 
problems,  the  use  of  crossover  exclusively  (along  with  a  suitable  population  size)  was 
shown  to  be  effective  (Deb  and  Agrawal,  1999:20).  This  is  not  surprising  since  crossover 
is  the  key  search  operator  and  implicitly  manipulates  the  best  substrings  or  building 
blocks  to  create  Pareto  optimal  solutions  (this  is  discussed  further  in  the  section  on  the 
Schema  Theorem).  This  does  not  imply  that  mutation  is  unimportant.  Mutation  is  used 
to  uncover  building  blocks  from  which  the  crossover  operator  may  direct  the  search  away 
(premature  convergence).  Too  large  a  rate  may  destroy  the  information  carried  by 
building  blocks.  Depending  on  how  much  pressure  the  researcher  wants  to  apply  to 
Pareto  front  distribution,  suggested  rates  for  mutation  range  between  0.001  and  0.1 
(Practical  Guide,  no  date).  De  Jong’s  work  with  GAs  on  problems  with  discontinuities, 
high  dimensionality,  noise,  and  multimodality  suggests  that  settings  of  n  =  50,  pc  -  0.60, 
and  pm  =  0.001  would  give  adequate  results  in  most  cases  (Mathematical  Optimization, 
no  date).  A  commonly  held  opinion  regarding  parameter  settings  is  that  “although  there 
is  no  unique  combination  guaranteeing  good  performance,  choosing  wisely  may  well 
result  in  more  effective  and  efficient  implementations”  (Van  Veldhuizen,  1999:2-18). 
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Quantitative  methods  for  determining  the  population  size  have  been  derived  in  the 
literature  but  Deb  and  Agrawal  have  recognized  that  what  is  needed  is  a  “good  yet  ready- 
to-use  population  sizing  estimate  for  generic  problems”  (1999:21).  For  crossover  based 
GAs,  Goldberg,  Deb,  and  Clark  derive  an  estimate  for  the  minimum  population  size,  Ns, 
needed  to  trigger  correct  building  block  processing  (Deb  and  Agrawal,  1999:1 1-12): 

2 

Ns=2ck ^8£-  (2.9) 

d 

where  c  is  the  tail  of  the  Gaussian  distribution  relating  to  the  permissible  error  rate  a,  k, 
is  the  number  of  competing  schemata,  and  crM2  /  d2  is  the  inverse  of  the  signal-to-noise  in 
the  underlying  problem. 

Both  De  Jong  and  Hessner  and  Manner  suggest  quantitative  methods  for 
determining  the  mutation  rate,  suggesting  that  the  rate  is  inversely  proportional  to  the 
population  size.  The  Hessner  and  Manner  formulation  is 


where  n  is  the  population  size  and  /  is  the  length  of  the  chromosome  (Practical  Guide,  no 
date). 

The  Schema  Theorem.  Chromosomal  representation  allows  the  manipulation  of 
information  about  the  search  space  and  the  transfer  of  that  information  to  other 
chromosomes.  This  information  is  carried  in  the  substrings  of  the  chromosome.  Thus, 
each  substring  represents  a  subspace  solution  (Osyczka  and  Kundu,  1995:95).  Substrings 
are  grouped  based  on  similarity  at  certain  string  positions,  called  schema,  and  are 
represented  on  a  template  of  0’s,  1  ’s,  and  *’s  (in  a  binary  representation)  (Dorigo  and 
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Maniezzo,  1993:8).  The  is  a  wildcard  symbol  that  represents  both  0  and  1 .  Thus  a 
schema  S\  =  (10**)  represents  stings  with  a  1  in  the  first  position  and  a  0  in  the  second 
position. 

The  length  of  a  schema,  8(S),  is  defined  as  the  distance  between  the  two  most 
distant  symbols  in  the  schema  that  are  not  wildcards.  The  order  of  a  schema,  o(S),  is 
defined  as  the  number  of  wildcards  subtracted  from  the  number  of  symbols  in  the  schema 
(Dorigo  and  Maniezzo,  1993:8).  So  S(S\)  is  2  -  1  =  1  and  o(S\)  is  4-2  =  2.  The  fitness 
of  a  schema  is  the  average  fitness  of  all  strings  that  match  the  schema  (Osyczka  and 
Kundu,  1995:95). 

Since  a  schema  is  a  grouping  of  similar  strings,  it  represents  a  region  in  the  search 
space.  For  the  objective  function  in  equation  (2.6),  the  schema  S\  represents  strings  with 
x  and  y  values  varying  from  8  to  1 1  with  function  values  varying  from  128  to  242.  A 
schema  S2  =  (00**)  would  result  in  function  values  varying  from  0  to  18.  Since  the 
objective  is  to  maximize,  strings  similar  to  S\  are  preferred  and  increase  in  proportion 
over  those  like  S2.  This  is  given  by  Holland’s  Schema  Theorem,  which  formalizes  the 
expected  number  m  of  schemata  h  within  a  population  at  generation  t: 


m{h,t  +  \)>m(h,t) 


m 


f  L 


8{h)  n, 

■  Pc-jff-  Pn,°(h) 


(2.11) 


where  flh)  is  the  average  fitness  of  all  strings  similar  to  h  within  the  population, /is  the 
average  fitness  of  the  population,  and  the  rest  as  defined  previously  (Goldberg,  et  al., 
1989:5).  It  can  be  seen  from  equation  (2.1 1)  that  short,  low-order,  above-average  fitness 
schemata,  or  building  blocks,  are  desirable  if  a  schema  is  to  grow  in  subsequent 
generations.  Building  blocks,  according  to  Goldberg,  “will  increase  in  number  with 
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exponential  speed”  (Dorigo  and  Maniezzo,  1993:9-10).  That  a  GA  accomplishes  this 
implicitly  through  the  selective  pressure  fostered  by  representation  schemes  and  genetic 
operators  is  postulated  by  Goldberg  in  what  is  known  as  the  Building  Block  Hypothesis 
(Deb,  1999a:14,  Van  Veldhuizen,  1999:4-3). 

Constraint  Handling.  Most  real  world  problems  are  going  to  be  constrained  in 
some  way  (time,  money,  space,  bandwidth,  etc.).  In  constrained  problems,  complexities 
arise  in  GAs  due  to  how  the  genetic  operators  direct  the  search.  It  is  very  likely  that  a 
small  change  to  a  feasible  solution  will  lead  to  an  infeasible  one.  (Ruiz-Andino,  et  al., 
2000:353). 

One  approach  to  dealing  with  constraints  is  to  modify  the  solution  representation 
itself  so  as  not  to  allow  the  creation  of  infeasible  solutions.  Repair  algorithms  or 
decoders  are  special  operators  that  avoid  the  construction  of  illegal  solutions.  They  may 
work  reasonably  well  but  are  highly  problem  specific  and  may  be  computationally 
intensive  to  run.  Additionally,  they  may  work  against  the  inherent  search  properties  of  the 
GA  and  may  be  difficult  to  implement  (Ruiz-Andino,  et  al.,  2000:353). 

Penalty  functions  that  reduce  the  fitness  of  infeasible  solutions  are  more  popular 
(Kundu,  1995:96,  Deb,  1999a:  14)  but  also  problem  specific.  Penalty  functions  maybe 
linear,  quadratic,  logarithmic,  etc.  functions  of  the  deviation  of  the  constraints  and/or  the 
number  of  violated  constraints.  Although  successfully  used  by  many  researchers,  the 
performance  of  GAs  will  depend  upon  the  choice  of  constraint  parameter  values  used. 

To  prevent  the  emphasis  of  a  particular  constraint  and  thereby  restrict  the  search, 
different  penalty  parameters  should  be  used  with  different  objective  functions  (Deb, 
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1998:7,  Gray,  et  al.,  1997).  Deb  describes  a  parameterless  penalty  function  that  is  used 
with  a  size-2  tournament  selection  operator  (Deb,  1999a:  15): 

Given  a  single  objective  function,  fix),  and  the  maximization  inequality 

gj(x)>  0,  y=l,2,...,  J  (2.12) 

the  fitness,  F(x),  of  any  solution  is  defined  as  follows: 


F(x)  = 


/(*)> 

J 

/max+Z^/W’ 

M 


ifgy(x)>0,  V /'gJ, 

otherwise 


(2.13) 


where  fmax  is  the  maximum  function  value  of  all  feasible  solutions  in  the 
population. 

In  a  tournament  between  an  infeasible  solution  and  feasible  solution,  it  can  be 
seen  from  equation  (2.13)  that  the  feasible  solution  always  has  a  better  fitness  than  the 
infeasible  ones.  If  both  solutions  are  feasible,  their  assessment  is  based  on  their 
respective  objective  function  values.  If  both  solutions  are  infeasible,  then  the  assessment 
is  made  based  on  the  amount  of  the  constraint  violations.  No  penalty  parameter  need  be 
used  since  pairwise  comparison  of  the  infeasible  solutions  does  not  depend  on  their  exact 
fitness  values  (Deb,  1999a:  15). 

Multiobjective  GAs.  The  basic  operation  of  the  single  objective  GA  in  Figure  3 
must  be  enhanced  to  evaluate  solutions  with  multiple  fitnesses  (objectives).  Researchers 
have  responded  with  a  number  of  ways  to  judge  the  overall  fitness  of  the  solutions.  Van 
Veldhuizen’s  recent  research  served  to  classify  known  multiobjective  EAs  on  the  basis  of 
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the  role  of  the  decision  maker  in  the  process  (1999:A-1).  In  a  priori  techniques,  the 
decision  maker  makes  his  or  her  preferences  known  at  the  beginning  of  the  process, 
resulting  in  a  single  compromise  solution,  e.g.  lexicographic  (ordering),  linear  and  non¬ 
linear  combination,  or  goal  programming. 

In  a  posteriori ,  the  Pareto  optimal  set  is  generated  for  the  decision  maker  who 
then  makes  his  or  her  preferences  known  by  selecting  a  solution  from  the  set.  This 
technique  is  notable  in  that  the  resultant  solution  set  is  independent  of  the  decision 
maker’s  preference  and,  assuming  no  change  in  the  problem  environment,  new  solution 
sets  would  not  need  to  be  generated  for  different  decision  makers. 

Progressive  techniques  allow  the  decision  maker  to  interact  and  provide 
preference  information  during  the  process.  These  techniques  require  a  high  degree  of 
participation  from  the  decision  maker  and  generally  make  use  of  both  a  priori  and  a 
posteriori  techniques. 

A  large  number  of  methods  forjudging  overall  fitness  use  an  objective- 
aggregation  approach  and  fall  in  the  category  of  a  priori  techniques.  The  different  fitness 
values  are  weighted  and  summed  according  to  the  decision  maker’s  preference  for  them. 
However,  this  is  very  subjective  and  is  difficult  to  do  accurately,  especially  when  the 
interplay  between  non-commensurate  objectives  is  not  well  understood  (Chipperfield,  et 
al.,  2000,  Shaw,  1998).  The  search  space  is  inextricably  linked  to  the  weightings,  thus  a 
single  inaccurate  weight  may  cause  a  GA  converge  to  an  undesirable  front. 

The  predominant  approach  to  solving  MOEAs  is  to  use  the  concept  of  Pareto 
dominance,  as  defined  in  Appendix  A,  in  the  selection  operator  (Deb,  1999:4).  Van 
Veldhuizen  notes  that  “the  sheer  number  of  Pareto  sampling  approaches  indicates  many 
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researchers  see  merit  in  the  basic  methodology”  (1999:3-10).  Pareto  dominance  allows 
all  nondominated  solutions  to  have  the  same  preference,  resulting  in  a  set  of 
nondominated  solutions  for  which  the  population-based  EA  is  particularly  well  suited  to 
handle.  Pareto  dominance  approaches  produce  as  their  end  product  nondominated  sets  of 
solutions  and  so  are  well  suited  for  use  in  the  a  posteriori  mode.  The  following 
algorithms  are  among  the  most  often  cited  and  copied  contemporary  MOEAs  that  use 
Pareto  dominance: 

Multiobjective  GA  (MOGA)  (Fonseca  and  Fleming,  1993) 

Niched  Pareto  GA  (NPGA)  (Horn,  et  al.,  1993) 

Non-dominated  Sorting  GA  (NSGA)  (Srinivas  and  Deb,  1995) 

Strength  Pareto  EA  (Zitzler  and  Thiele,  1998) 

NSGA-II  (Deb,  et  al.,  2000) 


In  order  to  improve  the  explorative  and  exploitative  properties  of  their  respective 
algorithms,  researchers  have  used  more  complex  selection  operators,  such  as  ranking, 
sharing,  niching,  elitist,  and  domination  tournaments  (Zydallis,  et  al.,  1999:2). 

MOGA  and  NSGA  use  variations  on  Goldberg’s  nondominated  sorting  procedure. 
The  basic  operation  of  this  procedure  is  to  rank  solutions  in  nondominated  order  with  the 
best  solutions  being  the  least  dominated.  These  fittest  solutions  are  given  higher 
probabilities  of  producing  more  offspring  (Bentley  and  Wakefield,  1999:7).  The  MOGA 
checks  the  population  and  assigns  a  rank  of  1  to  all  nondominated  solutions.  Each  of  the 
other  solutions  is  ranked  based  on  the  number  of  solutions  that  dominate  it  (Deb, 
1999b:5).  Fitness  is  assigned  based  on  linear  or  exponential  interpolation  (Van 
Veldhuizen,  1999:3-22). 
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NSGA  also  checks  the  population  and  assigns  a  rank  of  1  to  all  nondominated 
solutions,  forming  what  it  calls  the  first  level  of  non-domination.  The  first  level  is 
removed  from  consideration  and  proceeds  on  the  rest  of  the  population  in  the  same  way, 
resulting  in  1  to  n  domination  levels,  n  being  the  population  size.  All  first  level  solutions 
receive  a  fitness  equal  to  the  population  size.  The  other  levels  receive  a  dummy  fitness 
that  is  smaller  than  the  smallest  shared  fitness  of  the  preceding  level. 

Fitness  sharing  was  suggested  by  Goldberg  allow  for  solutions  with  identical 
fitness  along  different  parts  of  the  front,  thereby  helping  the  population  to  be  distributed 
along  the  front  (Horn,  et  al.,  1993:4).  The  number  of  neighboring  solutions  along  the 
front,  referred  to  as  a  niche,  are  used  to  selectively  reduce  the  fitness  of  high  niche  count 
solutions,  thus  increasing  pressure  toward  a  uniform  Pareto  front  distribution  (Horn,  et 
ah,  1993:8).  A  sharing  parameter,  oshare,  is  a  defined  maximum  distance  within  which 
any  solution  constitutes  as  belonging  to  a  neighborhood.  Solutions  within  <Jshare  of  each 
other  reduce  each  others  fitness. 

To  perform  domination  ranking,  NPGA  uses  domination  tournaments  of  size  two 
(Horn,  et  ah,  1993:  6).  The  tournament  procedure  selects  two  solutions  at  random  and 
each  of  them  competes  against  a  comparison  set  of  solutions,  tdom,  that  are  also  selected  at 
random  from  the  population.  When  one  solution  is  dominated  and  the  other  is  not,  the 
latter  is  selected.  When  competing  solutions  are  either  both  dominated  or  both 
nondominated,  sharing  determines  the  winner.  NPGA  implements  sharing  in  a  different 
manner.  Rather  than  reducing  the  fitness  of  high  niche  count  solutions,  the  winner  is 
declared  based  on  the  solution  with  the  smallest  niche  count  (Horn,  et  al.,  1993:9).  While 
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MOGA,  NSGA,  and  NPGA  all  require  explicit  values  for  crshare,  NPGA  also  requires  the 
same  for  tdom- 

Elitist  selection  ensures  that  the  best  solutions  are  retained  in  the  next  generation 
(Van  Veldhuizen,  1999:A-25).  SPEA  uses  an  elitist  selection  with  nondomination  (Deb, 
1999b:6).  The  algorithm  maintains  a  secondary  population  that  is  the  current  Pareto 
optimal  set.  This  population  is  combined  with  the  current  population  and  nondominated 
comparisons  are  performed  on  the  whole.  Nondominated  solutions  are  assigned  a  fitness 
based  on  the  number  of  solutions  they  dominate.  Preference  is  given  to 

1 .  nondominated  solutions  that  dominate  more  solutions  in  the  combined 
population,  and 

2.  dominated  solutions  that  are  dominated  by  more  solutions  in  the  combined 
population. 

The  preference  rules  are  meant  to  check  premature  convergence  by  preventing  large 
numbers  of  good  solutions  from  being  carried  over  from  one  generation  to  the  next  (Deb, 
1998:26).  However,  rather  than  having  equal  preference  for  all  nondominated  solutions, 
SPEA  is  biased  in  favor  of  nondominated  solutions  that  dominate  more  solutions  than 
others  (Van  Veldhuizen,  1999:3-20). 

There  is  another  approach  that  steps  further  away  from  the  traditional  GA:  the 
Messy  GA  (mGA).  The  mGA  abandons  fixed  length  bit  strings  and  so-called  neat 
operator,  crossover,  in  favor  of  variable-length  based  representations,  and  a  gene 
reordering  operator  called  cut-n-splice.  The  pseudo  code  for  Goldberg’s,  et  al.  mGA  is 
presented  in  Figure  5. 
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Start  mGA 

//  loop  for  user  defined  number  of  eras 
while  not  done  do 

// perform  Phase  1:  Partially  Enumerative  Initialization 
//  evaluate  individual  fitness  for  entire  population 
//  start  Phase  2:  Primordial  Phase 

//loop  for  user  defined  number  of  generations 
while  not  done  do 

//  increment  the  generation  counter 
g  =  g+  l; 

//  perform  Tournament  Thresholding  Selection 
//  test  for  appropriate  number  of  elapsed 
generations,  g* 
if  g  =  =  g*; 

reduce  population  size; 

//  reset  g 
g  =  0; 

end  if 

end  do 

//  end  Primordial  Phase 

II  start  Phase  3:  Juxtapositional  Phase 

II  loop  for  specified  number  of generations 
while  not  done  do 

//  increment  the  generation  counter 
t  =  t  +  1; 

//perform  cut-n-splice 

//  evaluate  individual  fitness  for  entire  population 
// perform  Tournament  Thresholding  Selection 

end  do 

//  end  Juxtapositional  Phase 
II  update  competitive  template 

end  do 
end  mGA 

Figure  5.  Pseudo  Code  for  a  Messy  Genetic  Algorithm  (Van  Veldhuizen,  1999:4-5) 


A  GA  has  difficulty  when  the  genes  are  not  ordered  properly.  According  to  the 
Schema  Theorem,  these  longer  length  schemata  have  a  higher  probability  of  being 
destroyed  by  crossover  and  mutation.  This  linkage  problem  leads  to  what  is  known  as 
deception,  where  poorly  ordered  schemata  lead  the  GA  away  from  the  global  optimum. 
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This  is  illustrated  in  the  following  example  by  Goldberg,  et  al.,  (1989:508).  Given  that 
(00****)  and  (****00)  are  highly  fit  schemata  of  an  optimal  point 
( 1  1  1  1  1  1  1  ),  and  that  the  schema  (00***00)is  much  less  fit  than  the  building 
block  (  1  1  *  *  *  1  1  ),  the  GA  will,  with  high  probability,  destroy  the  longer  length 
building  block  and  converge  to  a  less  than  optimal  point. 

Nature  allows  individuals  to  carry  redundant  information  such  as  multiple  copies 
of  genes  and  paired  chromosomes.  Messy  genetic  algorithms  copy  this  by  allowing 
redundant  or  even  contradictory  genes  (Goldberg,  et  al.,  1989:501).  To  allow  the 
reordering  of  genes,  each  gene  is  a  pair  of  integers  that  represents  the  name  and  value  of 
the  gene,  respectively.  For  example,  in  messy  representation,  two  strings  P\  and  P2  are 

/>  =  ( (3,1)  (1,0)  (3,0)  (2,1)  (1,0) )  (2.14) 

and 

P2=(  (4,1)  (2,0)  (3,0)  (2,1))  (2.15) 


Both  P\  and  P2  are  valid  despite  under-specification  by  P\  in  bits  four  and  five,  and  over¬ 
specification  by  P\  in  bit  3  and  by  P2  in  bit  2. 

The  traditional  crossover  operator  is  replaced  by  the  cut-n-splice  operator.  The 
name  of  this  operator  is  indicative  of  what  it  does  to  bit  strings.  The  position  of  cuts  can 
be  chosen  independently  for  both  parents.  After  cutting,  partial  strings  are  spliced  in  a 
random  order.  To  illustrate,  cut-n-splice  is  performed  on  Pi  and  P2  in  Figure  6. 

0\  (3,1)  (1,0)  (2,1) 

P,  (3,1)  (1,0)  (3,0)  (2,1)  (1,0) 

t  02  (3,1)  (1,0)  (4,1)  (2,0)  (3,0) 

03  (4,1)  (2,0)  (3,0)  (3,0)  (2,1)  (1,0) 


P2  (4,1)  (2,0)  (3,0)  (2,1) 

f  04  (2,1)  (3,0)  (2,1)  (1,0) 

Figure  6.  Example  of  Cut-n-Splice  Operation  (Hoffman,  1997) 
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Evaluation  of  variable  lengths  strings  is  problematic  since  under-  and  over- 
specified  strings  must  have  their  lengths  changed  to  fit  with  the  objective  function. 
Goldberg,  et  al.,  settled  on  a  simple  first-come,  first-served  process  to  handle  over¬ 
specification  (1989:501).  Since  this  method  does  not  rely  on  bitwise  fitness  for  its 
choice,  it  is  not  biased  to  toward  deceptive  schemata.  Goldberg,  et  al.,  successfully 
handle  under-specification  through  their  use  of  competitive  templates  that  fill  in  the 
unspecified  bits  in  an  under-specified  string  (1989:521).  A  competitive  template  is 
initialized  randomly  and  used  in  the  first  era.  Thereafter,  the  best  solution  in  the  current 
era  is  used  as  the  competitive  template  for  the  next  era,  and  so  on.  A  competitive 
template  that  is  itself  a  locally  optimal  solution  to  a  problem  “accentuates  salient  building 
blocks”  by  ensuring  that  their  fitness  is  better  than  that  of  the  template  (Goldberg  et  al., 


1990:417). 


An  mGA  proceeds  in  two  phases.  Prior  to  the  first  phase,  the  population  is 
initialized  so  that  it  completely  enumerates  building  blocks  of  a  given  length.  This 
process  is  referred  to  by  Goldberg,  et  al.  as  Partially  Enumerative  Initialization 
(1990:505).  The  population  sized  is  determined  based  on  the  highest  order  k  deceptive 
string  expected  in  the  problem: 


n  =  2 


( i\ 


\kJ 


(2.16) 


where  /  is  the  length  of  the  chromosome  (Goldberg,  et  al.,  1989:505). 

In  the  primordial  phase,  tournament  selection  is  performed  on  successive 
populations  “to  create  an  enriched  population  of  building  blocks  whose  combination  will 
create  optimal  or  very  near  optimal  strings”  (Goldberg,  et  al.,  1989:505).  To  avoid  an 
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apples  and  oranges  comparison  of  substrings  that  do  not  refer  to  the  same  subfunction, 
the  mGA  only  compares  substrings  that  are  similar  to  each  other  to  an  extent  that  is 
defined  by  a  threshold  number  of  genes  in  common.  The  threshold  parameter  is  defined 
by  Goldberg,  et  al.,  to  be 

0  =  ^1  (2.17) 

where  /  is  the  chromosome  length  and  l\  and  h  are  the  respective  substring  lengths 
(1989:426).  The  number  of  substrings  to  check  against  6  is  defined  by  a  shuffle  number, 
nSh,  equal  to  the  chromosome  length  (Goldberg,  et  al.,  1989:427). 

As  the  primordial  phase  proceeds,  the  population  size  is  reduced  roughly  in  half 
by  selection  at  specified  intervals  since  only  the  better  building  blocks  need  to  be 
maintained.  This  phase  represents  a  key  difference  between  a  mGA,  which  explicitly  and 
directly  manipulates  building  blocks,  and  other  EAs  which  settle  for  implicitly 
manipulating  building  blocks. 

At  the  conclusion  of  the  primordial  phase,  the  juxtapositional  phase  proceeds  as  a 
traditional  GA  would  on  a  fixed  population  size,  albeit  using  cut-n-splice  instead  of 
crossover  to  lengthen  the  substrings.  Thresholding  is  also  used  in  the  selection  operator. 
In  essence,  an  mGA  tries  to  gather  information  on  building  block  relationships  first,  then 
searches  for  better  solutions  (Kargupta,  no  date). 

In  his  PhD  dissertation,  David  Van  Veldhuizen  extended  the  mGA  to  handle 
MOPs  and  developed  the  Multi-Objective  Messy  Genetic  Algorithm  (MOMGA)  (1999). 
There  is  a  competitive  template  for  each  objective  that  is  at  first  randomly  initialized  and 
then  updated  with  the  previous  era’s  best  individual  for  that  objective  (Zydallis,  et  al., 
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2001 :4).  Like  Horn,  et  al.,  Van  Veldhuizen  augmented  the  tournament  selection  operator 
with  a  niching  strategy  to  increase  domination  pressure  (1999:4-14).  As  with  NPGA,  the 
MOMGA  requires  explicit  values  for  cxsiwre  and  tdom  to  control  domination  pressure.  The 
MOMGA  uses  Fonseca’s  suggested  method  to  determine  <jsiwre'- 

k  k 

N  =  -& - - - —  (2.18) 

CJ  share 

where  N  is  the  number  of  individuals  in  the  population,  A,  is  the  difference  between  the 
maximum  and  minimum  objective  values  in  dimension  i,  and  k  is  the  number  of  distinct 
MOP  objectives  (Van  Veldhuizen,  1999:6-11). 

The  MOMGA  also  maintains  and  updates  a  list  of  known  Pareto  optimal  solutions 
P known  with  Pareto  optimal  solutions  from  current  generation  P current  (Van  Veldhuizen, 
1999:4-17).  Since  dominance  determination  is  at  worst  an  n2  algorithm,  n  being  the  list 
cardinality,  it  is  done  on  Pknown  at  the  termination  of  the  program  to  prevent  the  MOMGA 
by  being  bogged  down. 

As  shown  by  equation  (2.16),  the  initial  population  grows  exponentially  as  the 
building  block  size  k  is  increased,  creating  a  computational  bottleneck,  i.e.  0(lk). 

Zydallis,  et  al.  reduce  this  bottleneck  using  a  probabilistic  approach  to  initialize  the 
population  (2001 :5).  Probabilistically  Complete  Initialization  (PCI)  creates  a  controlled 
number  of  building  blocks  of  size  k.  Building  Block  Filtering  (BBF),  which  replaces  the 
Primordial  phase,  alternately  reduces  string  lengths  by  randomly  deleting  bits  from  the 
strings  and  performs  selection  on  the  strings.  This  continues  according  to  a  user  specified 
schedule  of  alternations  until  the  strings  are  of  length  k.  The  Juxtapositional  phase 
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proceeds  as  before.  This  approach  probabilistically  ensures  that  all  of  the  best  building 
blocks  are  in  the  initial  population  and  results  in  initial  population  growth  on  order  of  the 
initial  string  length — 0(1)  (Goldberg,  et  al.,  1993:7). 

Using  these  ideas,  Zydallis,  et  al.  modified  original  MOMGA,  creating  a 
multiobjective  fast  messy  GA  (MOMGA-II).  Their  research  shows  the  MOMGA-II  to  be 
more  efficient  than  the  MOMGA  while  reutilizing  much  of  the  same  code.  The 
MOMGA-II  was  also  applied  to  the  same  test  suite  as  the  original  MOMGA,  achieving 
similar  results  but  with  fewer  juxtapositional  generations  (2001:10). 

Summary 

To  answer  the  first  research  question,  we  began  by  defining  what  a  multiple 
objective  programming  problem  is,  what  it  means  to  be  globally  optimal,  the  concept  of 
Pareto  dominance,  and  introduced  classical  approaches  to  solving  MOPs.  We  then 
presented  modem  methods  for  solving  MOPs  in  terms  of  two  approaches:  local  search  in 
objective  space  and  population-based. 

The  next  chapter  answers  the  second  and  third  research  questions  by  presenting  a 
multi  objective  model  formulation  and  solution  methodology  for  the  research  problem. 
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III.  Methodology 


Introduction 

The  previous  chapter  was  directed  at  the  first  research  question,  which  asks  for  a 
review  of  MOP  methodologies.  The  intent  of  this  chapter  is  to  answer  the  second  and 
third  research  questions: 

•  What  are  the  forms  of  the  decision  and  objective  spaces? 

•  How  is  the  selected  MOP  methodology  evaluated? 

This  chapter  also  defines  the  experimental  methodology  that  is  used  to  answer  the  fourth 
question:  evaluate  the  MOP  approach  used  to  solve  the  research  problem. 

This  chapter  begins  with  the  research  problem  model  formulation  that  is  based  on 
the  ALP  Pilot  Problem  (Swartz,  1999).  The  next  section  presents  the  MOP  formulation 
and  describes  the  construction  of  the  research  problem’s  decision  variables,  objective 
functions,  and  constraints.  This  is  followed  by  a  description  of  the  target  MOP  used  in 
the  model.  Next,  the  motivation  for  selection  of  a  specific  MOP  methodology  and 
objectives  for  its  evaluation  are  discussed.  The  last  two  sections  present  evaluation 
metrics  and  the  solution  methodology  used  for  this  research. 

Model  Formulation 

This  research  problem  is  modeled  on  the  ALP  Pilot  Problem  presented  by  Stephen 
Swartz  (1999).  The  reader  is  referred  to  Appendix  B  for  the  background  information 
used  to  construct  the  model. 
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The  MOP  Processing  Model  is  depicted  in  Figure  7.  The  centerpiece  of  the 
model  is  the  MOP  Tool.  The  MOP  Tool  is  to  be  an  application  of  a  MOP  solution 
methodology  selected  from  the  literature.  The  inputs  to  the  MOP  tool  allow  for  the  MOP 
formulation  discussed  in  the  next  section. 


[Specified  Resource  Levejst>  <CMRR  Task  Suitability  Matrix 

MOP 
TOOL 


I  Task  Preferencei^^>  <TjMRR  Weight  Cost  Matrix  I 


Figure  7.  MOP  Processing  Model 


The  output  requirement  of  the  model  is  to  present  to  the  war  planner  a  Pareto 
optimal  set  of  force  mixes  from  which  to  select  the  desired  force  mix.  Using  some 
decision  making  methodology,  the  planner  can  choose  the  desired  Mission  Ready 
Resource  (MRR)  set  from  each  Pareto  optimal  set  associated  with  an  inflection  point  on 
the  task  preference  vector  (as  shown  in  Figure  12  of  Appendix  B).  This  piecewise 
solution  represents  the  decision  maker’s  preferred  MRR  sets  for  a  given  combat 
capability. 
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These  assumptions  imply  that  a  solution  set  produced  at  a  given  resource  level  is 
a  subset  of  solution  sets  at  greater  resource  levels.  This  must  be  reflected  in  the 
piecewise  solution  and  is  accomplished  by  proceeding  along  the  preference  vector  from 
the  origin  to  the  last  resource  level  or  vice  versa.  This  is  illustrated  by  the  three  force  mix 
sets  in  Figure  8.  For  purposes  of  illustration,  each  set  is  considered  nondominated.  There 
are  nine  possible  force  mix  threads  that  progress  from  the  lowest  to  highest  specified 
resource  levels:  {1,  2,  5},  {1,  2,  6},  {1,  2,  7},  {1,  3,  5},  {1,  3,  6},  {1,  3,  7},  {1,  4,  5}, 

{1, 4, 6},  and  {1,  4,  7}.  However,  it  is  seen  by  inspection  that  only  three  of  these  meet 
the  subset  criteria:  {1,  2,  7},  {1,  3,  6},  and  {1,  3,  7}. 

“Thread  of  Progression” 


Resource  Level  =  60 

Force  Mix  5  Force  Mix  6  Force  Mix  7 


Fa 

fb 

30 

30 

Fa 

Fb 

15 

30 

Fa 

F„ 

13 

47 

Force  Mix  2 


F, 

fb 

16 

11 

Resource  Level  =  27 

Force  Mix  3 


bsh 

D 

Resource  Level  =  12 

Force  Mix  1 


Fa 

FB 

6 

6 

Force  Mix  4 


■3H 

Fe 

Q 

21 

Figure  8.  Thread  of  Progression 


It  is  tempting  to  handle  this  requirement  using  a  series  of  constraints.  To 
illustrate,  we  select  a  progression  direction  that  begins  with  Resource  Level  12.  We  then 
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choose  the  upper  bound  for  each  MRR  type  at  that  resource  level  to  be  the  starting  point, 
or  lower  bound,  for  the  next  resource  level.  This  leaves  {1,2}  and  {1,  3}  as  feasible 
threads.  Now  the  upper  bound  for  each  MRR  type  at  this  level  is  the  starting  point  for  the 
next  level.  This  leaves  only  {1, 2,  7}  and  {1,  3,  7}.  Although  {1,  3,  6}  meets  the  subset 
requirement,  it  is  deemed  infeasible.  Using  the  lower  bound  as  the  starting  point  for  the 
next  level  allows  the  infeasible  threads  {1,  2,  5}  and  {1,  2,  6}.  The  same  kind  of  problem 
exists  when  starting  from  the  highest  resource  level  and  progressing  downward. 

It  is  my  opinion  that  the  best  way  to  handle  the  construction  of  force  mix  threads 
is  through  a  post-processing  algorithm.  The  algorithm  operates  on  level-wise 
nondominated  sets  of  force  mixes  and,  starting  at  the  lowest  (or  highest)  resource  level, 
constructs  threads  iteratively,  taking  into  consideration  all  possible  feasible  threads.  Due 
to  time  constraints,  construction  of  this  algorithm  will  not  be  undertaken  in  this  research. 

MOP  Formulation 

Given  m  tasks  and  n  MRR  types,  the  solution  set  is  an  m  x  n  matrix.  A  matrix 
element  is  a  decision  variable,  xy,  that  represents  the  number  of  MRRs  of  type  j  allocated 
to  tasks  of  type  i.  Assuming  that  each  daily  task  is  satisfied  by  exactly  one  MRR,  and 
that  no  interactions  exist  between  differing  MRR  types,  then  the  suitability,  S,  for  all 
MRRs  is  defined  by 

n  m 

s=Z  2a*,  (3-0 

7=1  1=1 

where  a,j  is  the  suitability  of  MRR  j  for  Task  i  and  x,j  is  the  number  of  MRRs  j  allocated 
to  task  type  i. 
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The  requirement  that  all  tasks  i  =  1, n  must  be  satisfied  at  a  particular  resource 
level  (RL)  k  is: 

n 

RLtaskk  ■  =  xi  j  (3 .2) 

7=1 

Since  the  desired  capability  for  a  task  is  set  by  the  decision  maker  and  defined  to  be 
static,  the  left-hand  side  of  equation  (3.2)  is  an  equality  constraint. 

The  requirement  that  all  MRR  types  j  =  1,. .  .,n  do  not  exceed  their  available 
number  at  a  particular  resource  level  k  is 

m 

RLmrrJ,k^Hxij,k  (3-3) 

1=1 

In  this  problem,  the  decision  variables  are  allowed  to  take  on  any  non-negative  integer 
value  so  long  as  they  do  not  exceed  the  specified  resource  level.  Therefore,  the  left- 
hand  side  of  equation  (3.3)  is  an  inequality  constraint. 

The  maximum  number  of  sorties  per  day  for  a  particular  asset,  A,  is  given  by  its 
turn  rate,  t.  For  a  quantity  d  of  asset  A,  the  total  turn  rate  is 

TTRA  =  (d A)(turn  rate A)  (3-4) 

Given  that  A  has  P  configurations  corresponding  to  P  MRR  types,  the  upper  bound  for 
any  combination  of  the  P  MRR  types  is 

P  m 

rra,>X2>,,r  (3.5) 

r=\  i= 1 

For  example,  let  the  number  of  tasks  be  one,  and  let  P=  {1,  3, 4}  be  the  set  of  MRR 
types  that  correspond  to  asset  A.  If  the  number  of  A  is  one  and  the  turn  rate  of  A  is  two, 
then  following  decision  variable  values  are  possible: 


46 


xu 

Xl,3 

Xlt4 

1 

1 

0 

1 

0 

1 

0 

1 

1 

2 

0 

0 

0 

2 

0 

0 

0 

2 

1 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

Mathematically,  this  table  is  represented  as 

zix,, 

r~\  i= 1 

*1,+*13+Xl4<2 

It  is  difficult  to  determine  what  the  actual  logistical  footprint  is  for  a  given  asset 
set.  At  the  very  least,  it  is  clear  that  for  each  additional  asset  deployed,  there  is  a 
corresponding  increase  in  lift  cost  for  additional  resources,  e.g.  fuel,  munitions,  etc. 
Research  conducted  by  Matt  Goddard  suggests  that,  for  F-16s,  the  relationship  of  asset 
quantity  to  lift  resource  consumption  is  linear  (2001).  Assuming  that  lift  consumption  is 
linear  and  without  interaction,  the  weight  consumption,  W,  and  volume  consumption,  V, 
for  all  MRRs  are 

n  m 

I  V/''  V  (3-6) 

7=1  (=1 

and 

n  m 

V=Y.  V;vv. ,  (3.7) 

j= 1  1=1 
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where  and  Aj  are  the  weight  and  volume  consumed  by  a  single  MRR  j. 


The  form  of  the  suitability  maximizing  /  lift  minimizing  MOP  with  A  asset  types, 

n  MRR  types,  at  a  resource  level  k,  and  decision  variables  { xij,  Xjj, ... ,  xmn  }  is 

maximize: 

n  m 

S  =  H  tiauxu 

M  i= l 

minimize: 

n  m 

j=\  1=1 

n  m 

7=1  i= 1 

subject  to 

n 

yj  xi  j  =  RLtaskk  i  for  i  =  1  to  m 

i= i 

(3.8) 

m 

y  x •  j  <  RLmrrk  j  for  j  - 1  to  n 

i=i 

(3.9) 

Pa  m 

yyx^  <  TTRa  for  a  =  1  to  A,  and 

(3.10) 

r=l  i=l 


P„  =  number  MRR  types  for  a 
{ xi,i ,Xjj, ... , xmn }  are  non-negative  integers 


The  number  of  constraints  resulting  from  equation  (3.8)  is  equal  to  the  number  of 
tasks.  These  constraints  ensure  that  the  total  number  of  sorties  for  Task  i  is  exactly  the 
desired  capability  at  that  resource  level.  The  maximum  value  for  any  decision  variable  is 
found  by  using  equation  (3.8)  and  allocating  all  task  capability  to  one  MRR  type.  This 
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information  is  important  to  the  MOP  programmer  who  must  allocate  computer  memory 
to  hold  the  value  for  each  decision  variable. 

The  number  of  constraints  resulting  from  equation  (3.9)  is  equal  to  the  number  of 
MRR  types.  These  constraints  ensure  that  no  MRR  type  can  be  allocated  a  number  of 
sorties  that  exceeds  the  given  resource  level.  These  constraints  are  also  used  when  there 
are  restrictions  on  the  available  number  of  any  MRR  type,  e.g.  attrition  or  changes  in 
asset  turn  rate.  It  is  important  to  note  that  each  constraint  refers  to  a  single  MRR  type. 

Target  MOP 

The  MOP  Tool  inputs  in  Tables  2  through  5  provide  linkage  between  this  thesis 
and  concurrent  ALP  research,  and  create  a  search  space  large  enough  to  serve  as  a 
reasonable  test  of  the  MOP  Tool’s  utility  to  the  problem.  The  inputs  are  completely 
notional  but  not  entirely  arbitrary. 

The  number  of  tasks  in  Table  1  and  MRR  types  in  Table  5,  along  with  their  task 
suitabilities,  are  set  to  provide  proper  input  to  concurrent  research  (Filcek,  2001).  The 
suitabilities  reflect  notional  but  reasonable  values  that  clearly  differentiate  the  MRR 
types.  The  same  can  be  said  for  the  lift  consumption  values  in  Table  5. 

To  keep  the  number  of  task  capability  decisions  by  the  decision  maker  at  a 
reasonable  level,  five  resource  levels  in  Table  2  were  specified,  equating  to  15  separate 
task  preference  decisions.  These  preferences  are  reflected  in  Table  3.  When  the  ratios 
are  applied  to  their  respective  resource  level  values,  the  result  is  the  capability  matrix  in 
Table  4.  The  values  are  rounded  to  a  whole  number  so  that  the  sum  across  tasks  is  equal 
to  the  resource  level. 
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The  values  of  the  resource  levels  were  chosen  to  create  solution  spaces  of 
increasing  size.  Given  three  tasks  and  five  MRR  types  and  a  resource  level  of  300  sorties 
per  day,  the  worst  case  number  of  possible  force  mixes  is  approximately  9.72  x  1019  (by 
equation  (B.2)). 

For  the  target  MOP,  it  is  assumed  that  there  is  no  restriction  on  the  available 
number  of  any  MRR  type;  the  combined  lift  will  not  exceed  the  maximum  available  lift; 
no  attrition;  and  that  each  asset  has  one  associated  MRR  type,  i.e.  one  sortie  per  day. 
These  simplifying  assumptions  are  made  to  meet  research  time  constraints,  but  also  allow 
this  groundbreaking  research  an  opportunity  to  explore  the  basic  problem  complexity. 


Table  1.  Tasks 


INDEX 

NOMENCLATURE  j 

|  1 

Air-to-Air  (AA) 

2  .  1 

Air-to-Ground  (AG) 

1  .3  i 

Precision  Bombing  (PB) 

Table  2.  Resource  Levels  (RLs) 


INDEX 

RL  (sorties  per  day)  j  j 

1  16  J 

[  2 

32 

3 

75  | 

[  4 

150 

1  5  :l 

|  300  . j 

50 


Table  3.  Desired  Task  Capability  Ratios 


r 

PERCENT  TO  TASK  5 

INDEX 

AA 

AG 

PB 

1 

6° 

30 

10 

2 

30 

60 

10 

3 

25  60 

15 

!  4  ■ 

20 

50 

30 

5  : 

20 

30 

50 

Table  4.  Desired  Capability  Matrix 


j 

TASK  (sorties  per  day) 

i 

INDEX 

AA 

AG  ! 

PB 

DECISION  SPACE 
CARDINALITY 

1 

10 

5  i 

1 

630,630  | 

2 

10 

20 

2 

159,549,390 

3 

19 

45  S 

11  i 

«  2.56  xlO12 

4  1 

30 

75  ; 

45 

«  1.48x1 01 6  | 

1  5  j| 

60 

90  i 

'  . i 

!  150 

«  4.37  x  1019  | 

Table  5.  Task  Suitability  /  Lift  Consumption  Matrix 


INDEX 

MRR 

Type 

TASK  SUITABILITY 

LIFT  CONSUMPTION  j  j 

AA 

AG 

PB 

WEIGHT  ; 
(short 
tons j 

VOLUME 

(cubic 

feet) 

1 

Fa 

0.800 

o 

O 

o 

O 

o 

o 

20.2  j 

1650.0 

1  2 

Fb 

0.300 

0.800 

0.001 

28.5 

2475.0  J 

3 

Fc 

0.600 

0.600 

0.100 

35.7  1 

2887.5 

4 

B, 

0.001 

0.001 

0.800 

19.9 

O 

in 

o 

r- 

5 

b2 

0.001 

0.001 

0.400 

22.5 

2200.0  j 
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The  complete  MOP  formulation  is  as  follows: 

Decision  variables:  Number  of  MRR  j  assigned  to  Task  i  -  {xu  ...  Xjj} 
Maximize: 

S  -  0.8x,  ,  +  0.3x,  2  +  0.6x,  3  +  0.001x,  4  +  0.001x,  5  +  0.4x2]  +  0.8x22  +  0.6x23 
+0.00  lx,  4  +  0.001x25  +  0.001x31  +  0.001x32  +  0.1x33  +  0.8x34  +  0.4x35 

Minimize: 


(3.11) 


W  =  20.2(x, ,  +x2  I  +  x3l)  +  28.5(x,  2  +  x22  +  x32)  +  35.7(x13  +x23+x33) 
+19.9(xI4  +  x24  +  x3  4)  +  22.5(X|  5  +  x2  5  +  x35) 

V  =  1650(x, ,  +x21  +  x3I)  +  2475(x,  2  +  x22  +  x32)  +  2887.5(x13  +  x23  +  x33) 
+1705(x,  4  +  x24  +  x34)  +  2200(x,  5  +  x2  5  +  x35) 


(3.12) 


(3.13) 


Subject  to: 


x,  i,...,x35  >0 

(3.14) 

{x,i,...,x35}e/ 

(3.15) 

x,  i  +  X,  2  +  x,  3  +  x,  4  +  x,  5  =  RLtaskm] 

(3.16) 

X2,l  +X2,2  +  *2,3  +X2,4  +  *2,5  =  RLtClskm2 

(3.17) 

x3  ,  +  x3  2  +  x3  3  +  x3  4  +  x3  5  =  RLtaskmi 

(3.18) 

x,  ,+x2ii+x3  ,  <  RLmrrm  l 

(3.19) 

X,  2  +*2,2  +*32-  RLmrrm,2 

(3.20) 

xi  3  +  x2  3  +  *3  3  -  RLmrrm  3 

(3.21) 

*]  4  +*2,4  +*3,4  ^  RLmrrnA 

(3.22) 

x,  5  +  x25  +x3  5  <  RLmrrm  S 

(3.23) 

where  m  is  the  Resource  Level  index  for  the  current  problem. 
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Motivation  and  Objectives 

The  aim  of  this  research  is  to  select  and  evaluate  a  MOP  tool  that  outputs  to  a 
post-processor  an  acceptable  series  of  resource  level  constrained  force  mixes  in  a 
reasonable  amount  of  time,  and  is  scalable  to  more  realistic  problem  sizes.  In  their  basic 
forms,  simulated  annealing  and  genetic  algorithms  are  easy  to  understand  and  implement. 
Both  methods  can  show  convergence,  albeit  slowly,  when  good  solutions  now  are  better 
than  great  solutions  later. 

While  many  candidate  approaches  exist  in  the  literature,  an  evolutionary-based 
approach  is  ideally  suited  to  the  search  for  nondominated  sets  of  solutions,  particularly 
for  multidimensional  problem  domains  and  large  search  spaces.  The  balance  between 
exploration  pressure  and  exploitation  pressure  can  be  controlled  nearly  independently  in 
GA,  allowing  flexibility  in  design  (Deb,  1999a:l  1).  In  addition,  GA  populations  and 
operators  can  be  parallelized,  allowing  scaling  to  large  problems  by  using  multiple 
processors  to  reduce  overall  computational  time. 

The  choice  to  design  an  MOEA  for  the  problem  domain  or  to  modify  an  existing 
MOEA  to  incorporate  problem  domain  knowledge  is  a  practical  one.  All  methods  rely 
heavily  on  computing  skills  for  implementation,  and  the  time  allotted  for  this  research  is 
limited.  While  a  simple  GA  can  be  designed  for  the  problem  at  hand,  clever 
enhancements  are  needed  to  increase  the  efficiency  and  effectiveness  of  the  algorithm  an 
make  it  truly  useful.  The  contemporary  MOEAs  cited  in  the  literature  review  are  well 
designed,  based  on  modem  MOEA  theory,  and  have  found  use  in  many  applications  (Van 
Veldhuizen,  1999).  Van  Veldhuizen  identified  quantitative  metrics  to  objectively 
compare  four  such  MOEAs  in  terms  of  their  efficiency  and  effectiveness:  MOGA, 
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NSGA,  NPGA,  and  MOMGA.  The  reader  is  referred  to  Chapter  II  for  an  overview  of 
these  MOEAs.  Using  a  carefully  designed  test  suite  of  MOPs  that  emphasized  certain 
genotypical  and  phenotypical  characteristics  from  throughout  the  MOP  domain  (concave, 
convex,  connected,  non-connected,  scalable,  uniformity,  and  non-uniformity),  he 
concluded  that  the  MOMGA  compared  favorably  with  the  other  MOEAs  and  surpassed 
the  effectiveness  of  several  of  them  (1999:  5-5,  5-7,  7-23,  Zydallis,  et  al.,  2001 :1, 14). 
These  results  validated  the  use  of  explicit  building  block  MOEAs  in  the  MOP  domain.  A 
similar  experiment  by  Zydallis,  et  ah,  pitted  the  MOMGA-II  against  the  same  MOEAs 
and  concluded  that  the  MOMGA-II  achieves  results  similar  to  that  of  the  MOMGA,  but 
in  a  more  efficient  manner  (2001 : 14). 

The  algorithm  selected  as  the  MOP  tool  is  the  MOMGA-II.  Van  Veldhuizen 
comments  that  “although  no  guarantor  of  continued  success,  any  search  algorithm  giving 
effective  and  efficient  results  over  the  test  suite  might  be  easily  modified  to  target 
specific  problems”  (1999:5-7).  MOMGA-II  experimental  results  on  pedagogical 
problems  are  encouraging  and  its  architecture  is  readily  modifiable  for  the  test  problem’s 
high  dimensionality  and  heterogeneous  objectives.  The  Pareto  dominance  routine  is 
simpler  to  implement  when  the  objectives  are  either  all  minimization  or  all  maximization. 
Therefore,  the  target  MOP’s  first  objective  is  to  minimize  the  multiplicative  inverse  of 
equation  (3.11). 

MOMGA-II  use  of  memory  is  illustrated  in  Figure  9.  Memory  usage  increases 
linearly  with  string  length,  the  slope  of  which  is  dependent  upon  the  population  size.  The 
maximum  problem  length  for  the  target  problem  is  120  bits  (from  Table  4),  requiring 
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only  1 .36  kB  per  population  member.  A  suitably  equipped  PC  can  accommodate  much 
larger  problems. 


String  Length  vs.  Memory  Req’d 


Population  Size 


String  Length  (bits) 


2048 


1024 

512 


Figure  9.  MOMGA-II  String  Length  vs.  Memory  Required 


The  first  objective  is  to  compare  P,me  and  PF,rue  from  the  enumerative  approach  to 
Pknown  and  PFknown  from  the  MOMGA-II.  Since  enormous  computational  resources  are 
required  to  enumerate  Ptrue  and  PF,me  beyond  index  1  from  Table  4,  all  absolute 
comparisons  can  only  be  made  at  a  single  specified  resource  level.  To  explore  the 
performance  of  the  algorithm  on  the  target  MOP,  an  analysis  of  selected  MOMGA-II 
operating  parameters  can  determine  the  individual  impact  of  those  parameters  on  the 
algorithm’s  effectiveness. 

The  No  Free  Lunch  theorem  tells  us  that  on  average,  all  search  algorithms 
perform  equally  well  or  equally  bad  over  all  problems  (Wolpert  and  Macready,  1996:2). 
Therefore,  applying  a  search  algorithm  without  regard  for  problem  domain  knowledge 
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can  lead  to  performance  no  better  than  a  random  search.  Two  aspects  from  the  research 
problem  domain  to  consider  are  interrelationships  between  decision  variables,  decision 
variable  string  length,  and  application  of  side  constraints. 

Each  decision  variable  is  a  member  of  two  groups:  Task  Number  and  MRR  Type. 
The  strongest  relationship  is  the  Task  Number,  as  evidenced  by  equations  (B.l)  and 
(B.2).  The  strict  equality  of  the  decision  maker’s  task  preferences  means  that  for  one 
decision  variable  to  increase,  one  or  more  others  from  the  same  Task  Number  group  must 
decrease.  This  is  also  true  for  MRR  types,  but  to  a  far  lower  extent.  The  best  genotypic 
representation  for  related  decision  variables  would  be  one  that  fosters  this  relationship. 
The  Building  Block  Hypothesis  suggests  a  representation  that  encodes  them  in  close 
proximity  to  one  another.  This  idea  can  be  used  in  the  primordial  phase  of  the  MOMGA- 
II  by  using  building  blocks  of  size  equal  to  the  bitstring  representation  of  the  decision 
variables. 

The  traditional  binary  encoding  of  decision  variables  is  used  for  the  target  MOP. 
Since  the  decision  space  of  target  MOP  is  defined  to  be  in  7,  the  MOMGA-II  can  be 
designed  to  use  either  a  standard  bit  string  length  for  each  decision  variable  or  bit  string 
lengths  that  depend  on  the  size  of  each  decision  variable.  While  the  latter  economizes  on 
the  memory  required  for  each  vector  solution,  the  destruction  of  infomiation  by 
recombination  may  be  biased  toward  the  larger  representations  (assuming  that  the 
random  number  generator  output  is  normally  distributed).  Therefore,  the  representation 
length  for  each  decision  variable  is  defined  to  be  that  for  the  largest  decision  variable. 

The  algorithm  designer  typically  has  two  choices  when  using  side  constraints  to 
ensure  a  feasible  region.  If  the  constraints  are  applied  within  the  algorithmic  process  in 
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order  to  incorporate  some  problem  domain  knowledge,  solutions  that  are  found  to  be 
infeasible  are  penalized  or  thrown  out,  diminishing  their  impact  on  the  search.  This 
works  against  a  GA,  which  uses  the  fitness  landscape  information  carried  by  populations 
to  direct  the  search,  and  risks  premature  convergence.  If  the  constraints  are  applied  a 
posteriori,  information  from  the  entire  fitness  landscape  is  potentially  available  to  the 
algorithm.  The  resultant  unconstrained  nondominated  front  is  then  culled  of  infeasible 
members.  The  comparatively  larger  search  space  may  require  a  greater  number  of  fitness 
evaluations  to  converge  to  the  Pareto  front.  It  may  also  be  that  the  unconstrained  PFlrue 
does  not  contain  any  point  of  the  constrained  PF,rue. 

Since  the  best  approach  is  unknown,  and  given  limited  time  to  complete  this 
research,  the  second  objective  is  to  demonstrate  the  ability  of  the  applied  algorithm  to 
produce  level-wise  nondominated  force  mix  sets  as  defined  by  the  model  formulation. 
Since  the  constraints  in  equations  (3.8)  and  (3.9)  are  not  handled  explicitly  by  the 
algorithm,  the  decision  variable  domains  must  be  specified  explicitly  using  the  bitstring 
representing  each  decision  variable.  This  version  of  the  MOMGA-II  will  use  a  standard 
length  binary  representation  for  all  decision  variables. 

Preliminary  runs  of  the  MOMGA-II  on  the  target  MOP  using  index  1  of  Table  4 
and  the  parameters  listed  in  Table  7  generated  no  feasible  solutions.  Referring  to  Figure 
10,  the  objective  space  of  the  target  MOP  at  the  lowest  resource  level  defines  a  point 
mass  feasible  field.  This  is  a  particularly  difficult  problem  for  any  optimization  method. 
When  an  algorithm  finds  itself  at  any  feasible  point,  it  is  surrounded  by  infeasible  points 
of  varying  density.  There  is  no  guarantee  that  the  next  move  off  of  a  feasible  point  will 
result  in  another  feasible  point. 
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Figure  10.  Three  Views  of  the  Target  MOP  Tri-objective  Space  for  Resource  Level  1 

For  Resource  Level  1,  the  probability  of  selecting  a  feasible  point  is  given  by  the 
quotient  of  the  decision  space  cardinality  and  the  binary  value  of  the  chromosome  length 
(assuming  a  one-to-one  mapping  of  the  decision  space  to  the  objective  space): 

^^«5.47xl(T13  (3.24) 

Since  the  building  blocks  found  by  the  algorithm  are  the  result  of  converging  to 
the  infeasible  front,  the  use  of  initially  random  templates  was  discarded  in  favor  of  user 
selected  templates  that  are  corrected  to  be  feasible.  The  templates  are  defined  prior  the 


58 


start  of  the  run,  are  different  for  each  objective,  and  are  not  updated  with  the  best 
individual  after  each  era.  Preliminary  results  indicate  that  this  approach  effectively 
stimulated  the  algorithm’s  search  in  feasible  space. 


Performance  Measures 

If  PFtnie  is  known,  the  Final  Generational  Distance  metric  can  be  used  to 
characterize  how  “far”  PFknown  is  from  PFtme : 


(3.25) 


where  n  is  the  number  of  vectors  in  PFknown,  p  =  2,  and  d,  is  the  Euclidean  distance 
between  each  vector  and  the  closest  element  of  PFtrue  (Van  Veldhuizen,  1999:6-15). 

An  MOEA  adds  elements  to  PFknown  over  a  number  of  generations.  The  number 
of  nondominated  vectors  that  are  added  can  vary  depending  on  how  much  of  the 
objective  space  that  the  algorithm  is  allowed  to  explore.  The  Overall  Nondominated 
Vector  Generation  (ONVG)  metric  can  be  used  to  measure  how  “good”  an  MOEA  is  at 
generating  desired  solutions  given  a  fixed  fraction  of  search  space  to  explore  (Van 
Veldhuizen,  1999:6-18): 

ONVG^\PFk\  (3.26) 


Alone,  this  metric  says  nothing  about  the  quality  of  PFknown-  Even  if  an  MOEA’s  ONVG 
is  equivalent  to  the  number  of  vectors  in  PFtnie,  it  may  be  that  PF,rue  does  not  contain  any 
element  of  PFknown ■  The  Error  Ratio  metric  takes  this  into  consideration  in  order  to 
characterize  how  well  PFknown  converges  to  PF,rue : 
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(3.27) 


i- 

n 

where  n  is  the  number  of  vectors  in  PF known  and 

_  JO  if  vector  i,  i  =  (l,...,  n)e  PFmie , 

[l  otherwise 

(Van  Veldhuizen,  1999:6-14). 

Experimental  Design 

Computational  Environment.  The  MOMGA-II  is  written  in  ANSI  C  and  compiled 

using  the  Sun  Workshop  Compiler  version  C  4.2.  It  can  be  executed  in  any  UNIX  / 
LINUX  environment.  For  this  research,  the  execution  platform  is  a  Sun  Ultra  10 
workstation  equipped  with  a  450  MHz  processor,  1  GB  RAM,  and  Solaris  2.8. 

The  target  MOP  is  a  discrete  mapping  from  the  integers  to  M3 .  On  a  computer, 
real  number  representation  is  dependent  on  machine  specific  resolution.  Therefore, 
determination  of  PF,nie  is  machine  dependent.  However,  for  the  three  objectives — 
suitability,  weight,  and  volume — accuracy  on  the  order  of  10'6  and  beyond  is  not  an  issue. 

Exhaustive  deterministic  enumeration  is  needed  to  find  Ptrue  and  PFtrUe ,  limited  by 
machine  specific  resolution,  of  course.  A  program  that  performs  this  task  was  written  in 
ANSI  C  and  compiled  using  Microsoft  Visual  C++  6.0.  The  source  code  for  this  program 
is  presented  in  Appendix  C.  The  program  code  was  not  optimized  and  uses  only  a  single 
list  to  maintain  P known  and  PFknown.  The  program  was  executed  on  a  Dell  Precision  610 
equipped  with  a  Pentium  II  500MHz  processor,  2GB  RAM,  and  using  Windows  2000 
Professional.  It  was  observed  that  this  platform  initially  processes  on  average  40,000 
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solutions  per  minute  but  that  this  rate  decreases  exponentially  with  each  additional 
solution.  Using  this  program,  complete  enumeration  of  scenarios  beyond  index  1  in 
Table  4  would  exceed  the  time  requirement  for  this  thesis  by  centuries.  Therefore, 
absolute  comparisons  are  only  made  between  the  enumerated  P,rue  /  PFlme  and  the 
MOMGA-II P known  /  PF known  for  index  1  in  Table  4.  The  enumerated  PFtrue  and  Ptrue  is  in 
Appendix  D. 

Experiments.  The  first  objective  is  to  compare,  in  terms  of  absolute  performance, 
the  output  of  the  MOMGA-II  to  the  baseline  exhaustive  enumeration  solution  in  order  to 
tune  the  algorithm  to  the  MOP  formulation.  There  is  no  definitive  answer  as  to  whether 
this  comparison  should  be  made  in  decision  or  objective  space.  This  decision  is  typically 
left  to  personal  preference.  For  this  research,  objective  space  comparisons  are  preferred 
because  of  the  high-dimensionality  of  the  decision  space.  For  statistical  comparison 
purposes,  all  results  will  be  taken  from  30  replications  of  the  MOMGA-II  using  a 
different  random  number  generator  seed  each  time.  The  randomQ  function  in  the 
random,  c  file  generates  a  single  random  number  between  0.0  and  1.0  using  the 
subtractive  method  described  in  (Knuth,  1981). 

The  second  objective  is  to  demonstrate  the  ability  of  the  implicitly  constraining 
MOMGA-II  to  produce  level-wise  nondominated  force  mix  sets  as  defined  in  the  section 
on  model  formulation.  The  third  objective  is  to  examine  how  execution  time  responds  to 
target  MOP  problem  size.  These  three  objectives  are  met  through  the  following 
experiments — absolute  performance  response  to  parameter  changes,  execution  timing, 
and  the  demonstration  of  level-wise  nondominated  force  mix  sets. 
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Absolute  Performance  Response  to  Parameter  Changes.  The  goal  of  this 
experiment  is  to  identify  how  MOMGA-II  effectiveness  changes  with  different  key 
parameter  settings.  There  are  no  studies  that  show  which  parameters  and  what  values  are 
key  to  good  performance  for  MOEAs  (Van  Veldhuizen,  1999:6-7).  While  a  complete 
parameter  analysis  is  warranted  in  this  case,  the  allotted  time  for  this  research  is  limited 
and  only  certain  primary  operators  and  values  can  be  investigated. 

Experimental  Parameters.  The  relationship  between  decision  variables 
may  affect  algorithm  performance.  Different  building  block  ranges  will  be  used  to  take 
advantage  of  this  possible  relationship.  In  the  MOMGA-II,  search  is  primarily 
accomplished  through  building  block  filtering,  splicing,  and  selection.  The  probability  of 
cutting  a  string  will  be  changed  to  investigate  its  affect  on  algorithm  effectiveness.  One 
of  these  settings  will  be  a  zero  probability  of  cutting,  placing  the  burden  of  search  on 
building  block  filtering,  splicing,  and  tournament  selection.  The  splice  operator  used  in 
the  juxtapositional  phase  is  the  primary  method  of  string  composition  in  the  MOMGA-II. 
Its  probability  will  be  reduced  to  reveal  its  affect  on  performance.  Finally,  initial 
population  size  affects  the  amount  of  fitness  landscape  information,  and  therefore 
building  blocks,  available  to  building  block  filtering.  PFtrue  will  be  compared  to  the 
results  of  the  MOMGA-II  by  individually  applying  each  of  the  parameter  settings  listed 
in  Table  6.  Table -6r-to  the  base  settings  listed  in  Table  7.  A  total  of  nine  alternative 
parameter  settings  are  used,  each  with  30  replications.  The  metrics  G,  ONVG,  and  E  are 
used  to  quantitatively  compare  the  alternatives. 
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Table  6.  Experimental  Parameter  Settings 


|  PARAMETER 

SETTINGS 

BB  size 

{2,4,8} 

|  Pcul 

{0,0.2,  .2} 

|  P  splice 

{  1.0,  0.8,  0.6  } 

init.  pop  size 

{  600,915,  1200} 

Table  7.  MOMGA-II  Parameter  Settings 


PARAMETER  i 

SETTING  j 

Pcul  i 

0-2  . | 

| 

P  splice  ! 

1.0  j 

. Pm  , . _ . _  j 

. . 0.0 . . 

tdom 

3  .  . . | 

&shtire  | 

Equation  (2.18) 

eras 

. 4 . J 

Initial  population  size 

915  . | 

termination  1 

average  string  length  «  problem  size  | 

Execution  Timing.  This  analysis  does  not  take  into  account  the  additional 
run  time  needed  to  obtain  an  acceptable  level  of  convergence  and  looks  only  at  how 
execution  timing  is  affected  by  increasing  the  problem  size.  For  each  index  in  Table  4, 
and  using  the  base  parameter  settings  in  Table  7,  a  time  hack  will  be  taken  from  program 
start  to  the  end  of  the  juxtapositional  phase  in  era  4.  Thirty  replications  are  used  for  each 
index  to  provide  a  good  statistical  sample. 

Demonstration  of  Level-wise  Nondominated  Force  Mix  Sets.  The  goal  of 
this  experiment  is  to  demonstrate  the  ability  of  the  MOMGA-II  to  produce  level-wise 
nondominated  force  mix  sets  in  the  absence  of  explicit  constraint  handling  methods. 
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Since  Ptnie  and  PFtme  are  unknown  for  scenarios  beyond  index  1,  no  statistical 
comparisons  are  made.  Of  particular  interest  is  non-dominated  set  cardinality  after 
applying  the  constraints. 

MOMGA-II  Parameter  Settings.  Unless  otherwise  stated,  each  replication  of  the 
MOMGA-II  is  performed  using  the  settings  listed  in  Table  7.  Except  for  the  termination 
rule,  the  listed  parameter  settings  are  the  “default”  settings  used  in  previous  research 
(Van  Veldhuizen,  1999,  Zydallis,  et  al.  2001b).  The  termination  rule  was  selected  to 
allow  runs  with  larger  building  block  sizes  to  be  less  dependent  on  template  fitness. 

Summary 

This  chapter  addressed  the  second  research  question  by  presenting  a  model  to 
describe  the  inputs  and  outputs  to  a  MOP  tool  and  formulating  the  research  problem 
MOP  in  terms  of  its  decision  variables,  objective  functions,  and  constraints.  The  target 
MOP  was  constructed  using  inputs  that  provide  linkage  between  this  thesis  and 
concurrent  ALP  research,  and  also  create  a  search  space  of  sufficient  size  with  which  to 
evaluate  the  MOP  tool. 

After  discussing  the  selection  of  the  MOMGA-II  as  the  MOP  tool,  the 
experimental  objectives,  metrics,  and  key  MOMGA-II  parameter  settings  were  presented 
in  order  to  address  the  third  research  question:  how  to  evaluate  the  selected  MOP 
methodology.  The  chapter  concluded  with  a  discussion  of  the  computational 
environment  and  the  three  experiments  designed  to  meet  experimental  objectives.  The 
next  chapter  addresses  the  fourth  research  question  by  analyzing  MOMGA-II 
effectiveness  and  efficiency  when  applied  to  the  target  MOP. 


64 


IV.  Results 


Introduction 

The  previous  chapter  outlined  the  experimental  methodology  used  to  apply  the 
MOMGA-II  in  three  comparative  experiments.  This  chapter  focuses  on  the  last  research 
question,  which  addresses  MOMGA-II  solution  quality  and  efficiency  as  it  pertains  to  the 
research  problem  of  optimizing  MRR  suitability  and  lift  cost.  The  sections  that  follow 
present  the  analysis  of  the  experimental  data  and  report  the  results. 

First,  the  results  from  parametric  testing  form  the  basis  for  an  absolute 
comparison  of  the  MOMGA-II  to  PFtnie  for  a  single  resource  level.  Next,  execution 
timing  results  are  used  comment  on  the  time  complexity  of  the  MOMGA-II  as  it  relates 
to  problem  scale.  Finally,  the  results  from  level-wise  runs  are  used  to  demonstrate  the 
ability  of  the  implicit  constraint  handling  MOMGA-II  to  produce  non-dominated  sets  of 
force  mixes  corresponding  to  various  sortie  resource  levels. 

Statistical  Analysis 

Absolute  Performance  Comparison.  The  sample  data  for  Final  Generational 
Distance,  Figure  1 3  in  Appendix  F,  suggests  that  the  population  distributions  are  not 
normal,  requiring  non-parametric  methods  for  statistical  comparisons.  The  technique 
employed,  the  Kruskal-Wallis  H- test,  assumes  that  the  samples  are  random  and 
independent,  at  least  five  measurements  in  each  sample,  and  that  the  probability 
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distributions  from  which  the  samples  are  drawn  are  continuous  (McClave,  et  al., 
1998:892).  These  assumptions  are  satisfied  and  so  the  following  hypotheses  are  tested: 

Ho:  The  probability  distributions  of  the  parameter  groups  for  the  G  metric 

are  the  same. 

HA:  At  least  two  of  the  groups  are  different. 

The  H- test  in  Figure  13  of  Appendix  F  was  accomplished  using  JMPIN  4.0.2.  At 
the  0.1  significance  level,  the  observed  significance  level  of  0.15  indicates  that  there  is 
insufficient  evidence  to  reject  Ho. 

The  ONVG  sample  data  across  the  parameter  groups  are  only  borderline  normally 
distributed  and  so  the  Kruskal- Wallis  H- test  is  also  applied  here  (Figure  14  in  Appendix 
F).  The  ONVG  sample  data  are  cardinal  and  therefore  Poisson  distributed  (Reynolds, 
2001).  For  a  large  ONVG  range  (in  this  case  from  0  to  630,630)  the  assumption  of  a 
continuous  distribution  is  reasonable  (Reynolds,  2001).  The  following  hypotheses  are 
tested: 

Ho:  The  probability  distributions  of  the  parameter  groups  for  the  ONVG 

metric  are  the  same. 

Ha:  At  least  two  of  the  groups  are  different. 

At  the  0.1  significance  level,  the  observed  significance  level  of  0.29  indicates  that 
there  is  insufficient  evidence  to  reject  H0. 

These  results  show  that  using  implicit  constraint  handling,  MOMGA-II 
convergence  properties  and  non-dominated  set  cardinality  are  not  significantly  affected 
by  the  choice  of  selected  parameter  values  within  the  range  tested.  When  infeasible 
solutions  are  not  explicitly  dealt  with  by  the  algorithm,  it  is  relying  primarily  upon 
building  block  filtering  and  the  juxtapositional  generations  to  select  and  combine  good 
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building  blocks  into  feasible  solutions.  Target  MOP  pilot  runs  of  the  MOMGA-II 
without  any  modifications  to  the  algorithm  resulted  in  no  feasible  solutions  at  all.  By 
simply  specifying  a  feasible  template  for  each  objective  and  carrying  them  forward  every 
era,  non-dominated  set  cardinality  improved  to  a  mean  of  5  with  a  standard  deviation  of 
0.4. 

Referring  to  Figures  15  -  23  in  Appendix  G,  the  three-dimensional  plots  of  PF,rue 
and  PF known  for  each  parametric  alternative  allow  qualitative  assessments  of  MOMGA-II 
performance.  Although  the  MOMGA-II  never  found  any  solutions  in  PFtrue,  the  PFknoWn, 
nearly  all  of  the  non-dominated  sets  approximate  the  structure  of  PFtrue  and  some  of  the 
solutions  are  very  close  to  the  front.  That  this  occurred  without  the  incorporation  of  any 
explicit  constraint  handling  methods  is  indicative  of  the  algorithm’s  robustness. 

Execution  Timing  Analysis.  This  analysis  does  not  take  into  account  the 
additional  run  time  needed  to  obtain  acceptable  level  of  convergence,  looking  only  at 
how  execution  timing  is  affected  by  increasing  the  problem  size.  In  all  cases,  the 
MOMGA-II  completed  4  eras  in  under  25  seconds.  It  is  seen  from  Figure  1 1  that  within 
the  range  of  60  to  120  bits,  execution  time  increases  nearly  linearly. 

Without  extrapolating,  it  is  difficult  to  see  how  execution  time  responds  to  scaling 
up  the  test  MOP  to  a  real  world  size  problem  with  seven  basic  aerospace  missions  and  on 
the  order  of  1 00  or  more  MRR  types.  This  preliminary  result  suggests  that  on  the  tested 
platform,  MOMGA-II  execution  timing  may  scale  linearly  or  possibly  quadratically  with 
chromosome  length. 
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Demonstration  of  Level-wise  Nondominated  Force  Mix  Sets.  The  feasible 
solution  set  cardinalities  for  the  two  lowest  resource  levels  are  both  six.  The  MOMGA-II 
found  no  nondominated  feasible  solutions  at  any  of  the  other  resource  levels.  The 
fractions  of  feasible  points  to  total  points  corresponding  to  each  resource  level  differ  from 
each  other  by  approximately  an  order  of  magnitude: 


RL, 

RL2 

rl3 

rl4 

rl5 


6.30  x  105 

260 

1.60  x  108 

275 

2.56  x  1012 

o90 


1.48  x  10 

o105 


16 


4.37  x  10 
/-)120 


,19 


«  5.47  x  1 0"13 


«  4.22  x  10'15 


«  2.07  x  10"15 


«  3.65  x  10'16 
«  3.29  x  10‘17 
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The  fractions  of  feasible  space  are  not  grossly  smaller  at  progressively  larger 
resource  levels  and  may  not  completely  account  for  the  algorithm’s  inability  to  find 
feasible,  nondominated  solutions  above  resource  level  two.  With  random  population 
initialization,  the  algorithm  may  be  misled  right  off  the  bat  by  starting  with  infeasible 
members.  Nor  does  the  use  of  feasible  templates  guarantee  that  feasible  building  blocks 
will  be  found  during  building  block  filtering.  The  same  feasibility  problems  apply  to  the 
operations  in  the  juxtapositional  phase.  Since  the  templates  are  not  updated  with  the  best 
found  individual  for  each  objective,  this  version  of  the  MOMGA-II  operates  with  a 
handicap.  The  overall  result  is  that  the  algorithm  is  drawn  to  an  infeasible  front. 

A  recently  updated  version  of  the  MOMGA-II  incorporates  a  solution  repair 
function  that  iteratively  adjusts  bits  in  a  random  manner  until  a  solution  becomes  feasible, 
restricting  the  search  to  feasible  space.  The  target  problem  was  applied  to  the  algorithm, 
but  due  to  time  constraints,  proper  tuning  and  complete  analysis  was  not  performed. 
Preliminary  results  are  promising.  Using  the  same  base  case  parameters,  random  seeds, 
and  30  replications,  the  reported  metrics  in  Table  8  are  greatly  improved. 

Table  8.  Explicit  Constraint  Metrics 


MEAN  j 

STANDARD 

DEVIATION 

MEDIAN 

G 

5.2 

7.7 

0.1 

ONVG 

22.7  ; 

1.5 

23 

1  E 

49.8 

i.i 

47.9 

When  results  are  combined  over  4  replications,  the  MOMGA-II  finds  all  26  points  in 
PFtrue  for  Resource  Level  1 . 
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When  applied  to  Resource  Levels  1  through  5,  the  MOMGA-II  was  able  to 
produce  feasible,  equally  preferred  force  mixes  with  the  cardinalities  listed  in  Table  9. 
Three-dimensional  plots  of  PFknown  for  each  level  (Appendix  F)  suggest  that  a  solution 
basic  structure  is  maintained  over  all  problem  scales. 

Table  9.  Explicit  Constraint  MOMGA-II  ONVG  by  Resource  Level. 


RESOURCE 

LEVEL 

ONVG  i 

1 

26 

2  1 

70 

3  i 

131 

4  | 

215 

5 . j 

301  | 

Summary 

This  chapter  answers  research  question  four  by  presenting  experimental  analyses 
and  results  for  three  experiments  that  address  the  experimental  objectives:  absolute 
performance  response  to  parameter  changes,  execution  timing,  and  the  demonstration  of 
level-wise  nondominated  force  mix  sets.  Overall,  the  results  reveal  the  MOMGA-II’s 
robustness  and  linear  execution  time  (on  the  range  tested),  and  show  that  implicit 
constraint  handling  as  applied  to  the  MOMGA-II  does  not  go  far  enough  to  prevent 
convergence  to  an  infeasible  front.  Preliminary  results  of  an  explicit  constraint  handling 
version  of  the  MOMGA-II  show  greatly  improved  performance. 
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V.  Conclusion 


Introduction 

Chapter  I  began  with  the  motivation  of  this  research,  stating  that  the  Defense 
Advanced  research  Projects  Agency’s  Advanced  Logistics  Project  (ALP)  seeks  to  bring 
campaign  planning  into  the  21st  century  with  multiple,  real-time  deployment  plans  and 
rapid  replanning.  We  can  capitalize  further  on  what  ALP  brings  to  the  table  by  providing 
a  methodology  to  evaluate  multiple  plans  and  provide  the  warfighter  with  the  best 
available  package  with  which  to  do  the  job. 

As  a  front-end  to  ALP,  the  Mission-Resource  Value  Assessment  Tool  (M-R 
VAT)  intends  to  assess  competing  force  mixes  in  terms  of  their  intrinsic  task  capability 
and  the  campaign  specific  issues  affecting  their  effective  employment  in  the  theater.  The 
primary  goal  of  this  research  was  to  identify  force  mixes  that,  in  terms  of  their  intrinsic 
value,  represent  the  best  match  of  assets  to  tasks  with  the  smallest  deployment  footprint. 

To  accomplish  this  goal,  four  research  questions  answered: 

1 .  Which  methodologies  can  be  used  to  trade-off  intrinsic  value  and  deployment 

cost  (lift)  and  result  in  a  set  of  force  mixes  that  are  preferred  over  others? 

2.  What  are  the  forms  of  the  decision  and  objective  spaces? 

3.  How  is  the  selected  approach  evaluated? 

4.  Does  the  selected  approach  result  in  an  acceptable  solution  to  the  research 

problem  in  a  reasonable  amount  of  time? 

A  series  of  research  phases  was  used  to  answer  these  questions. 
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The  first  phase  was  a  literature  review  intended  to  provide  a  broad  overview  of 
the  multiobjective  optimization  problem  (MOP)  class  in  order  to  help  with  the  selection 
of  an  appropriate  methodology.  The  review  highlighted  concepts  of  MOP  formulation, 
competing  objectives,  global  versus  local  optimization,  constrained  optimization,  and 
Pareto  dominance  of  solutions.  Also  reviewed  were  classical  and  modem  methodologies, 
including  those  from  the  metaheuristic  class.  The  results  of  the  literature  review  served 
to  answer  research  question  1 . 

The  next  phase  of  this  research  incorporated  the  ALP  Pilot  Problem  (Swartz, 
1999)  in  the  definition  of  the  multiobjective  model  and  formulation.  The  Mission  Ready 
Resource ,  defined  as  a  building  block  of  capability,  suitability,  and  deployment  cost;  the 
Task  Preference  Vector,  and  the  definitions  of  acceptable  force  mixes  defined  the 
decision  space,  objective  space,  and  the  problem’s  constraints,  thus  answering  research 
question  2.  This  phase  also  defined  the  target  multi  objective  problem  used  to  evaluate 
the  MOP  model  and  formulation,  and  their  application  to  the  selected  MOP  methodology. 

An  important  goal  of  this  research  was  to  incorporate  as  much  problem  domain 
knowledge  as  possible  into  the  algorithmic  approach.  The  third  phase  used  the  problem 
domain  knowledge  from  the  model  and  MOP  formulation  to  select  an  appropriate  MOP 
methodology.  The  Multi  objective  Messy  Genetic  Algorithm  (MOMGA-II),  shown  to 
robustly  and  flexibly  handle  a  variety  of  difficult  pedagogical  problems,  was  selected  to 
produce  the  desired  nondominated  sets  of  solution,  thus  answering  research  question  3. 
The  algorithm  was  adapted  to  the  level-wise  output  requirement  of  the  MOP  model.  The 
MOMGA-II  also  incorporates  implicit  constraint  handling  in  order  to  investigate  that 
method’s  affect  on  the  effectiveness  of  the  algorithm. 


72 


Finally,  the  last  phase  evaluates  the  employed  methodology  based  on  three 
experimental  objectives.  The  first  objective  was  to  compare,  in  terms  of  absolute 
performance,  the  output  of  the  MOMGA-II  to  the  baseline  output  of  the  exhaustive 
enumeration  solution  in  order  to  help  tune  the  algorithm  to  the  MOP  formulation.  The 
second  objective  was  to  demonstrate  the  ability  of  the  implicitly  constraining  MOMGA- 
II  to  produce  level-wise  nondominated  force  mix  sets  as  defined  by  the  model 
formulation.  The  third  objective  was  to  examine  how  execution  time  responds  to  target 
MOP  problem  size.  The  answer  to  research  question  4  is  found  in  the  experimental 
results  showing  that  implicit  constraint  handling  as  applied  to  the  MOMGA-II  is  not  a 
viable  approach  to  producing  acceptable  sets  of  nondominated  force  mixes.  However, 
the  results  also  reveal  the  MOMGA-II’s  robustness  and  linear  execution  time  (on  the 
range  tested).  This,  along  with  the  greatly  improved  solution  quality  and  cardinality 
achieved  by  preliminary  runs  of  the  solution  repairing  MOMGA-II  support  the  viability 
of  this  approach  to  producing  well  balance  force  mixes. 

Conclusions 

This  research  shows  that  the  multiobjective  model  and  problem  formulations, 
along  the  with  test  problem,  are  constructive  approaches  to  investigating  the  problem  of 
force  mix  selection.  The  importance  of  this  research  is  in  the  illumination  of  problem 
complexity  for  so  abstract  a  problem.  Although  simplified  by  assumptions,  the  employed 
methodology  allows  us  to  gauge  problem  complexity  and  uncover  problem  domain 
knowledge  that  must  be  incorporated  into  any  solution  platform.  An  important  aspect  of 
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problem  complexity  was  revealed  by  mathematically  defining  the  cardinality  of  the 
constrained  multiobjective  research  problem  decision  space. 

A  program  allowing  for  deterministic  enumeration  and  Pareto  dominance 
checking  of  a  large  solution  space  was  developed  to  support  the  research  methodology. 
In  addition,  a  post-processing  program  was  developed  to  analyze  MOMGA-II  output  for 
solution  dominance,  feasibility,  and  uniqueness.  Both  programs  were  essential  in 
supporting  this  methodology  and  have  applications  beyond  it  as  part  of  a  generic  MOP 
tool  kit. 

While  implicit  constraint  handling  allows  the  MOMGA-II  complete  access  to 
fitness  landscape  information,  this  research  demonstrated  that  such  a  method  is 
ineffective  when  the  unconstrained  Pareto  front  does  not  contain  the  constrained  Pareto 
front.  Despite  this  handicap,  this  research  showed  the  adaptability  of  the  MOMGA-II  to 
the  problem  domain  and  viability  of  the  platform  for  larger  scale  problems.  Preliminary 
results  indicate  that,  without  parametric  tuning,  the  solution  repair  function  significantly 
contributes  to  algorithm  convergence  to  feasible,  optimal  force  mixes,  warranting  further 
investigation  as  an  efficient  and  effective  method  for  identifying  well  balanced  force 
mixes. 

Limitations 

A  number  of  simplifying  assumptions  were  made  for  this  research:  single  sortie, 
capability,  linear  suitability,  and  no  limit  on  available  lift  or  assets.  These  assumptions 
dull  the  real-world  applicability  of  the  devised  model  and  must  be  dealt  with  before 
employing  the  model  outside  of  the  research  realm. 
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The  execution  time  experiment  was  conducted  on  a  chromosome  length  range 
that,  until  now,  had  not  been  attempted.  The  reality  of  the  situation,  however,  is  that  the 
number  of  USAF  MRR  types  ranges  in  the  hundreds  and  those  of  our  allies  compound 
this  situation.  While  the  “approximately  linear”  result  of  the  experiment  is  positive, 
extrapolation  of  the  algorithm’s  time  complexity  to  a  real-world  scaled  problem  is  fraught 
with  peril. 

For  practical  reasons,  parametric  tests  of  the  algorithm  were  limited  in  terms  of 
range  and  number  of  parameters.  This  constitutes  only  a  coarse  tuning  and  is  probably 
not  the  best  combination  of  parameter  values  to  use.  As  revealed  by  the  literature  review 
on  genetic  algorithms,  parameter  settings  are  typically  problem  dependent  and  little 
conclusive  research  on  proper  parameter  settings  exists.  Parameter  value  selection  is 
more  art  than  science  at  this  point. 

A  design  consideration  of  the  M-R  VAT  front-end  is  that  it  should  run  on  a 
Microsoft  Windows  enabled  PC  so  that  it  may  be  widely  distributed  throughout  the 
planning  community.  At  this  time,  the  MOMGA-II  code  executes  on  Unix  and  Linux 
platforms.  Significant  effort  is  needed  to  port  the  code  to  run  on  a  Microsoft  platform. 

As  part  of  the  methodology,  problem  domain  knowledge  was  sought  out  and 
applied  to  the  MOP  tool  when  found.  There  may  be  elements  of  the  problem  that  have 
escaped  the  eye  of  the  researcher,  elements  which,  when  properly  employed  can  effect 
the  effectiveness  and  efficiency  of  the  search  algorithm. 

The  model  used  in  this  research  is  based  strictly  upon  the  intrinsic  value  of  force 
mixes.  The  intrinsically  scored  output  of  the  model  is  then  extrinsically  scored  to 
produce  a  ranking  of  force  mixes.  However,  force  mixes  that  have  composite  intrinsic 
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and  extrinsic  values  may  have  a  different  Pareto  front,  much  less  a  different  fitness 
landscape.  Extrinsic  knowledge  constitutes  problem  domain  knowledge  and  can  be 
incorporated  as  part  of  the  MOP  formulation.  Another  tie-breaking  methodology  would 
need  to  be  employed,  but  this  too  may  also  be  incorporated  as  problem  domain 
knowledge.  The  choice  of  where  to  stop  with  this  reasoning  may  best  depend  on  how 
much  problem  domain  flexibility  is  lost  or  gained  by  the  tool — a  highly  accurate 
assessment  means  nothing  if  the  rules  for  assessment  change  and  the  tool  is  unable  or  too 
complex  to  adapt. 

Recommendations 

This  research  was  limited  in  that  only  an  implicit  constraint  handling  method  was 
employed  in  the  MOMGA-II.  Any  complete  assessment  of  the  algorithm’s  utility  to  the 
research  problem  must  include  an  explicit  constraint  handling  version.  Preliminary 
results  from  the  MOMGA-II  using  a  solution  repair  function  showed  greatly  improved 
solution  quality  and  cardinality,  clearly  indicating  the  direction  of  future  efforts  to 
identify  well  balanced  force  mixes. 

When  Ptrue  or  PFm,e  is  unknown,  some  of  the  quantitative  metrics  employed  by 
Van  Veldhuizen  can  assess  convergence  to  PFtrue  (1999:8-5).  Relative  quantitative 
performance  reveals  differences  between  alternatives  but  offers  little  in  terms  of  solution 

quality.  The  best  assessment  is  an  absolute  comparison  of  a  MOP  tool’s  solution  output 

2 

to  P,rue  or  PFtrue.  Complete  enumeration  of  all  possible  solutions  is  limited,  as  it  is  an  n 
algorithm,  n  being  the  number  of  solutions  in  the  space.  It  is  recommended  to  completely 
enumerate  at  least  three  inflection  points  (sortie  levels)  on  the  Task  Preference  Vector  in 
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order  to  quantitatively  evaluate  a  selected  algorithm’s  solution  quality  and  move  beyond 
empirical  results.  The  IBM  SP  computers  at  both  the  Aeronautical  Systems  Center’s 
Major  Shared  Resource  Center  (ASC  MSRC)  and  the  U.S.  Army  Corps  of  Engineers 
Waterways  Experiment  Station’s  (CEWES)  MSRC  can  be  used  to  deterministically 
enumerate  all  possible  solutions  for  a  given  MOP  at  levels  that  are  well  beyond  the 
computational  capabilities  of  mere  desktop  machines. 

Future  Research 

For  the  MOP  formulated  in  this  research,  it  is  possible  that  the  structure  of  PFtrue 
is  not  sensitive  to  problem  scale.  The  advantage  in  this  case  is  that  the  search  can  be 
localized  to  the  region  of  the  objective  space  that  PFtrue  lies  for  any  number  of  decision 
variables  and  sortie  level.  As  suggested  in  the  recommendations  section,  it  would  be 
beneficial  to  completely  enumerate  several  points  on  the  Task  Preference  Vector  to 
explore  this  possibility. 

The  research  problem  can  also  be  applied  to  two  or  more  different  MOP  solving 
approaches,  all  of  which  can  be  compared  to  the  current  deterministic  approach  in  terms 
of  flexibility,  solution  quality,  ease  of  implementation,  and  scalability. 

In  terms  of  problem  scale,  chromosome  length  is  a  limiting  factor  to 
computational  efficiency.  Efficiency  can  be  gained  by  employing  variable  length  strings 
so  that  the  length  of  the  bit  string  representing  a  complete  decision  variable  is 
independent  of  the  others.  This  will  reduce  computer  memory  overhead  and  the  size  of 
the  search  space,  but  possibly  at  the  expense  of  bias  in  favor  of  operating  on  longer 
length  decision  variables. 
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One  way  to  handle  the  construction  of  force  mix  threads  is  through  a  post¬ 
processing  algorithm.  The  algorithm  operates  on  level-wise  nondominated  sets  of  force 
mixes  and,  starting  at  the  lowest  (or  highest)  resource  level,  constructs  threads  iteratively, 
taking  into  consideration  all  possible  feasible  threads.  However,  the  worst  case 
maximum  number  of  threads  to  construct  given  K  resource  levels  is 

K 

nkXnMX-"XnK=Ylnk  (5-1) 

*= 1 

where  n  is  the  number  of  MRRs  to  expend  at  a  given  resource  level.  Such  a  potentially 
large  number  of  alternative  force  mixes  would  overwhelm  a  decision  maker  and  must 
therefore  be  pared  to  a  reasonable  number.  A  way  to  do  this  that  considers  the  best 
feasible  threads  is  to  use  Filcek’s  extrinsic  scoring  model  to  rank  the  individual  force 
mixes  for  each  resource  level  (2001).  Rather  than  randomly  selecting  a  starting  point  and 
any  following  points,  construction  of  feasible  force  mix  threads  proceeds  using  the 
highest  ranked  force  mixes  first.  Not  all  of  the  possible  feasible  force  mix  threads  need 
be  constructed,  but  those  that  are  constructed  are  based  on  the  best  intrinsically  and 
extrinsically  evaluated  force  mixes. 

Finally,  the  simplifying  assumptions  used  to  initially  explore  the  research  problem 
may  almost  certainly  have  a  large  impact  upon  the  problem  domain  and  the  fitness 
landscape.  It  is  recommended  that  the  need  for  these  assumptions  be  evaluated  and, 
where  needed,  supplanted  by  more  realistic  modeling  information  in  order  to  enhance  its 
operational  applicability. 
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Summary 

The  multiobjective  optimization  problem  model  and  formulation  developed  for 
the  Mission-Resource  Value  Assessment  Tool  establish  the  fundamental  form  and 
complexity  of  force  mix  selection  defined  by  the  ALP  Pilot  Problem.  The  research 
methodology  and  its  results  are  the  first  steps  to  providing  rapid  force  mix  selection  based 
on  task  suitability,  sortie  capability,  and  the  amount  of  finite  lift  resources  consumed. 
Follow-up  research  is  well-positioned  to  spring  forward  from  this  point  and  develop  the 
multi  objective  force  mix  assessment  tool  capable  of  rapidly  providing  best  value  /  small 
footprint  alternative  force  mixes  corresponding  to  the  desires  of  the  warfighter. 


79 


Appendix  A:  Pareto  Concepts 


MOPs  present  a  set  of  solutions  from  a  trade-off  surface  between  objectives.  To 
be  better  understood,  it  is  necessary  to  define  key  Pareto  concepts:  Pareto  Dominance, 
Pareto  Optimality,  the  Pareto  Optimal  Set,  and  the  Pareto  Front.  Van  Veldhuizen’s 
definitions  are  noted  for  their  consistency  with  theory  and  are  the  ones  used  in  this  thesis 
(1999:2-3): 

Pareto  Dominance:  A  vector  u  =  (ux,...uk)  is  said  to  dominate  v  =  (v,,...,vt) 

(denoted  by  u  <  v )  if  and  only  if  u  is  partially  less  than  v ,  i.e., 

for  all  /  e  {1, . . . ,  k) ,  ui  <  v;  and  there  exists  an  i  e  {1, . . . ,  k }  such  that  ui  <  v; . 

Pareto  Optimality:  A  solution  x  e  Q  is  said  to  be  Pareto  optimal  with  respect  to 
Q  if  and  only  if  there  is  nor'eQ  for  which  v  =  F(x’)  =  (ft  (x'), . . . ,  fk  (V)) 

dominates  ii  =  F(x)  =  (fl(x),...,fk (x)) .  The  phrase  “Pareto  optimal”  is  taken  to 
mean  with  respect  to  the  entire  decision  variable  space  unless  otherwise  specified. 

Pareto  Optimal  Set:  For  a  given  MOP  F(x),  the  Pareto  optimal  set  (P")  is 
defined  as: 

P *  :={igD  such  that  there  is  no  xeQ  where  F(x')  -<  F(x)}  (A.l) 

% 

Pareto  Front:  For  a  given  MOP  F(x)  and  Pareto  optimal  set  P  ,  the  Pareto  front 
(PF*)  is  defined  as: 

PF*  :=  {u  =  F(x)  =  (/ (x), ...,fk (x))  |  x e  P' }  (A.2) 

P  is  the  set  of  all  solutions  whose  vectors  are  non-dominated  with  respect  to  all 
other  vectors  in  the  objective  space.  When  mapped  to  the  objective  space,  the  solutions 
in  P  form  the  Pareto  front,  PF. 
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Appendix  B:  Advanced  Logistics  Program  (ALP)  Pilot  Problem 

A  Mission  Ready  Resource  (MRR)  is  a  combination  of  an  asset  type  and  its 
resources,  e.g.  aircraft,  pilot,  fuel,  munitions,  support  equipment  and  personnel,  etc.,  that 
is  designed  to  have  a  certain  suitability  for  a  single  task.  MRR  suitability  to  a  given  task 
is  measured  per  sortie  on  an  absolute  scale  from  0  to  1 ,  with  0  indicating  no  suitability 
and  1  indicating  perfect  suitability  (Johnson,  2001).  A  combination  of  MRR  types  is 
defined  to  be  a  MRR  set  or  force  mix.  To  demonstrate,  assume  that  a  notional  aircraft  F, 
has  two  configurations.  Fa  and  F$,  which  constitutes  two  MRR  types.  Further,  if  the 
aircraft  could  be  prepared  and  flown  three  times  per  day,  then  it  would  represent  three 
MRRs  per  day.  These  three  MRR’s  could  either  be  all  FA  configurations,  or  all  FB 
configurations,  or  some  combination  of  the  two  configurations.  Finally,  assume  three 
tasks:  Suppression  of  Enemy  Air  Defenses  (SEAD),  Air-to-Air  (AA),  and  Combat  Air 
Support  (CAS).  Assuming  negligible  interaction  between  assets,  the  preference 
relationship  between  the  assets  and  tasks  is  illustrated  in  Table  8. 

Table  8.  Asset-Mission  Task  Preference  Matrix  (A-M  TPM)  1  (Swartz,  1999:1) 


The  critical  physical  dimensions  of  an  asset  are  its  pallet  weight  and  pallet 
volume.  A  single  MRR  requires  a  given  pallet  weight  and  volume  to  support  its 
deployment  and  use  for  a  given  sortie.  The  additional  complexity  of  translating  from 
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MRR  space  to  asset  space  is  beyond  the  scope  of  this  research,  hence  the  simplifying 
assumption  that  a  MRR  is  equivalent  to  a  single  asset  that  is  sampled  without 
replacement,  i.e.  flies  a  single  sortie  per  day. 

We  can  consider  a  MRR  as  a  building  block  of  1)  task  suitability,  2)  pallet  weight, 
and  3)  pallet  volume.  Task  suitability,  pallet  weight  and  pallet  volume  are  hereafter 
referred  to  simply  as  suitability,  weight,  and  volume,  respectively. 

A  realistic  assumption  is  that  the  decision  maker’s  task  preference  may  change 
depending  on  the  number  of  resources  at  his  /  her  disposal  and  the  phase  of  the  campaign. 
At  the  termination  of  a  campaign,  many  if  not  all  deployed  resources  are  expected  to 
redeploy  to  their  origins  or  some  other  location.  The  ALP  Pilot  Problem  is  concerned 
strictly  with  the  deployment  planning  and  execution  phase.  So  it  also  assumed  that 
resource  levels  do  not  decrease,  thus  modeling  a  build-up  of  resources  over  time;  and  that 
the  decision  maker  does  not  prefer  fewer  sorties  over  time.  For  instance,  given  a 
relatively  low  sortie  generation  level  at  the  beginning  of  a  campaign,  a  combatant 
commander  whose  goal  is  air  dominance  may  prefer  a  ratio  of  60  percent  AA,  30  percent 
SEAD,  and  only  10  percent  CAS.  Over  time,  the  sortie  generation  level  increases  and  the 
next  campaign  phase  may  emphasize  ground  attack.  This  is  reflected  in  the  combatant 
commander’s  task  preferences:  20  percent  AA,  30  percent  SEAD,  and  50  percent  CAS. 
The  decision  maker’s  task  preference  is  explicitly  indicated  by  the  number  of  sorties,  i.e. 
capability,  desired  for  each  task  at  a  given  resource  level  (sorties  available).  This  leads 
naturally  to  a  monotonic  task  preference  vector  with  components  that  sum  to  the  resource 
level.  The  points  at  which  the  ratio  of  desired  capability  changes  are  the  inflection 
points.  This  is  shown  in  Figure  12. 
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Vector  of  Mission  Preferences  over  Resources 
(Sorties  achievable  increases  along  vector) 


Each  MRR  Set: 

Collection  of  Mission-Ready  Resources 
satisfying  a  given  level/mix  of  sorties 
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Figure  12.  Matching  Resources  to  Tasks  (Johnson  and  Swartz,  2000b:  16) 


At  the  termination  of  a  conflict,  most  if  not  all  deployed  resources  are  expected  to 
redeploy  to  their  origins  or  some  other  location.  The  ALP  Pilot  Problem  is  concerned 
strictly  with  the  deployment  planning  phase.  Two  assumptions  are  made  in  the 
construction  of  the  preference  vector.  First,  it  is  assumed  that  resource  levels  do  not 
decrease,  thereby  modeling  a  build-up  of  resources  over  time.  Second,  the  preference 
vector  is  monotonic,  i.e.  the  decision  maker  does  not  prefer  fewer  missions  over  time. 

We  define  a  campaign  phase  as  the  time  period  where  the  commander’s  mission 
preferences  remain  constant.  If  the  commander  would  decide  that  fewer  missions  were 
required — in  a  possibly  different  ratio  of  task  types — then  a  new  campaign  phase  would 
begin  and  the  planning  process  described  here  would  be  repeated. 

A  decision  maker  can  make  his  /  her  task  capability  preferences  known  at  every 
feasible  resource  level,  thereby  minimizing  interpolation  error.  A  model  with  seven  tasks 
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and  a  maximum  resource  level  of  30  would  require  210  task-resource  allocation 
decisions,  unreasonably  overburdening  the  decision  maker.  Recognizing  this,  the 
Mission-Resource  Value  Assessment  Tool  (M-R  VAT)  model  uses  a  reasonable  number 
of  inflection  points  to  define  the  preference  vector  and  interpolates  preferences  for 
intermediate  points.  The  resulting  interpolation  error  is  assumed  to  result  in  a  negligible 
difference  between  the  preference  vector  and  the  decision  maker’s  true  task  preferences. 

The  decision  maker  has  a  finite  number  of  force  mix  choices  with  which  to  meet 
his  /  her  task  preferences.  The  number  of  force  mixes  for  a  single  task  is  given  by  the 
following  equation: 


(Jh  +m- 1)! 
nk  !(m  —  1) ! 


(B.l) 


where  nk  is  the  desired  sortie  level  for  Task  k,  and  m  is  the  number  of  MRR  types.  For  p 
tasks,  the  number  of  force  mixes  is 


nc,-,  <b-2> 

k-\ 

Using  the  example  in  Figure  12,  a  sortie  mix  of  25  AA,  20  CAS,  and  45  SEAD  results  in 
approximately  5.35xl013  unique  force  mixes  for  the  decision  maker  to  consider.  A  look 
at  the  number  of  aircraft  choices  and  configurations  available  in  today’s  Air  Force  makes 
it  clear  that  the  number  of  MRR  types  ranges  in  the  thousands  and  that  the  number  of 
possible  solutions  grows  unmanageably  large. 

Since  an  MRR  is  a  building  block  of  suitability,  weight,  and  volume,  the  choice  of 
the  best  MRR  set  depends  on  how  it  optimizes  the  three  criteria:  maximize  suitability, 
minimize  weight,  and  minimize  volume.  A  MRR  set  is  acceptable  if,  for  a  given 
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resource  level,  the  decision  maker  is  indifferent  to  the  tradeoff  between  the  three  criteria 
of  suitability,  weight,  and  volume.  Therefore,  the  formulation  of  acceptable  force  mixes 
is  a  multiobjective  problem  with  competing  objectives. 
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Appendix  C:  Source  Code  for  ENUMERATIONS 


/*  ENUMERATION. C  v2 . 4 

*  Author:  Dave  Wakefield 

*  email:  wakester@earthlink.net 

*  Date:  2/6/01 

*  Reference:  Author,  "IDENTIFICATION  OF  PREFERRED  OPERATIONAL  PLAN 
FORCE  MIXES  USING  A 

MULTIOBJECTIVE  METHODOLOGY  TO  OPTIMIZE  RESOURCE 
SUITABILITY  AND  LIFT  COST," 

Masters  Thesis.  Air  Force  Institute  of  Technology, 
Wright -Patterson  AFB ,  OH. 

*  2001. 

*  This  program  performs  a  complete  enumeration  of  a  decision  space 
defined  by  15  decision 

*  variables.  The  output  is  a  text  file  containing  in  column  format  the 
nondominated 

*  front  for  three  objectives  defined  in  the  evaluate  function,  and  the 
corresponding 

*  decision  variables.  Objective  1  is  maximized,  the  others  are 
minimized.  The  program 

*  takes  as  input  a  report  filename  and  whether  screen  output  of 
progress  is  desired.  A 

*  progress  file  named  'progress.txt'  is  created  and  updated  every  100th 
solution . 

*/ 


#include 

#include 

#include 

#include 

#include 


<stdio . h> 
<stdlib . h> 
ctime . h> 
<math. h> 
<conio . h> 


double  *Pcurrent ;  //  Archive  for  Pareto  front 

int  sizePcurrent ,*  //  Number  of  solutions  in  Pcurrent,  not  the 

index;  includes  solutions  flagged  for  removal 

double  fitness  [3]  =  {  0,  0,  0  };  //  initialize  fitness  vector 


void  evaluate (  double  f [] ,  int  il,  int  i2,  int  i3,  int  i4,  int  i5,  int 
i6 ,  int  i7 ,  int  i8,  int  i9,  int  ilO,  int  ill,  int  il2,  int  il3,  int  il4, 
int  il5  ) ; 

int  pareto (  int  n  ) ; 
double  factorial (  double  a  ) ; 
double  mod(  double  a,  int  b  ); 
void  clear  kb (  void  ); 


main  () 

{ 

int 

int 

int 

int 

int 


Taskl ;  //  DM  task  1  preference 
Task2;  //  DM  task  2  preference 
Task3 ;  //  DM  task  3  preference 

maxXll;  //  Upper  bound  for  Xll  (task  1,  MRR  type  1) 
maxX12;  //  Upper  bound  for  X12  (task  1,  MRR  type  2) 
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int 

maxX13 ; 

// 

Upper 

bound 

for 

X13 

(task 

i, 

MRR 

type 

3) 

int 

maxX14 ; 

// 

Upper 

bound 

for 

X14 

(task 

i, 

MRR 

type 

4) 

int 

maxXIS ; 

// 

Upper 

bound 

for 

X15 

(task 

i, 

MRR 

type 

5) 

int 

maxX21 ; 

// 

Upper 

bound 

for 

X21 

(task 

2, 

MRR 

type 

1) 

int 

maxX22 ; 

// 

Upper 

bound 

for 

X22 

(task 

2, 

MRR 

type 

2) 

int 

maxX23 ; 

// 

Upper 

bound 

for 

X23 

(task 

2, 

MRR 

type 

3) 

int 

maxX24 ; 

// 

Upper 

bound 

for 

X24 

(task 

2, 

MRR 

type 

4) 

int 

maxX25 ; 

// 

Upper 

bound 

for 

X25 

(task 

2, 

MRR 

type 

5) 

int 

maxX31; 

// 

Upper 

bound 

for 

X31 

(task 

3, 

MRR 

type 

1) 

int 

maxX32 ; 

// 

Upper 

bound 

for 

X32 

(task 

3, 

MRR 

type 

2) 

int 

maxX33 ; 

// 

Upper 

bound 

for 

X33 

(task 

3, 

MRR 

type 

3) 

int 

maxX34 ; 

// 

Upper 

bound 

for 

X34 

(task 

3, 

MRR 

type 

4) 

int 

maxX35 ; 

// 

Upper 

bound 

for 

X35 

(task 

3, 

MRR 

type 

5) 

int 

numTasks 

= 

3;  // 

Number 

of 

tasks 

int 

numMRRs 

5;  //  Number 

of  MRR  types 

int 

i;  //  used 

to  print  report 

int 

il,  i2. 

i3,  i4 , 

i5 ,  i6, 

i7 

18 , 

i9 ,  ilO, 

ill 

i!2 

/ 

il3 , 

il4,  il5  ; 

//  used  in 

for  statements 

int 

n;  //  Counts  solutions 

for  Pareto  dominance 

testing 

int 

z;  // 

expression  for 

switch 

statement 

int 

nondominated;  // 

1  if  solution 

is  nondominated,  0 

otherwise 

int 

t;  //  progress  report  solution 

increment 

time_t 

start,  finish,  time2t; 

// 

start 

,  finish, 

elapsed 

time 

// 

double 

aa,  bb, 

cc 

/ 

double 

duration; 

/ /  program  execution 

time 

double 

cardDS ; 

// 

Calculated  cardinality  of 

the  decision 

space 

double 

countDS 

= 

0;  // 

Counter 

for  number  of 

solutions 

evaluated 

double 

progress ; 

//  percentage 

of 

cardDS  evaluated 

FILE  *fpl ; 

FILE  *  f p2 ; 

char  filename [20] ; 

char  *progressf ile  =  "progress . txt " / 

char  output,  retry; 

/*  Print  program  description  to  screen  */ 

printf ("ENUMERATION. C  v2 . 4\nAuthor :  Dave  Wakef ield\nemail : 
wakester@earthlink . net\n" ) ,* 

printf ( "Date :  2/6/0l\nRef erence :  Author,  IDENTIFICATION  OF 
PREFERRED  OPERATIONAL  PLAN"); 

printf  ("FORCE  MIXES  USINGA  MULTIOBJECTIVE  METHODOLOGY  TO  OPTIMIZE 
RESOURCE  SUITABILITY  "); 

printf ("AND  LIFT  COST,  Masters  Thesis.  Air  Force  Institute  of 
Technology,  ") ; 

printf ( "Wright-Patterson  AFB ,  OH.  2001 . \n\nThis  program  performs 
a  complete  enumeration  of  ") ; 

printf ("a  decision  space  defined  by  15\ndecision  variables.  The 
output  is  a  text  file  "); 

printf ( "containing  in  column  format  the\nnondominated  front  for 
three  objectives  defined  in  "); 

printf  ("the  evaluate  function,  and\nthe  corresponding  decision 
variables.  Objective  1  is  maximized,  the  others  are  minimized.  The  "); 
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printf ( "program  takes  as  input  user-defined  sortie  levels  for 
Tasks  l\nthrough  3,  a  report  filename,  and  whether  screen  output  of 
progress  is  " ) ; 

printf ("  desired. A  progress  file  named  'progress.txt'  is  created 
and  updated  every  100th\nsolution . \n\n" ) ; 

/*  Get  user  defined  task  levels.  Allows  user  reinput  levels  as 
desired.  */ 

while  (  1  )  { 

puts (  "\nEnter  number  of  sorties  desired  for  Task  1"  ); 
scanf (  "%d",  &Taskl  ); 
clear__kb  ( )  ; 

puts (  "\nEnter  number  of  sorties  desired  for  Task  2"  )/ 
scanf (  "%dM,  &Task2  ); 
clear_kb ( ) ; 

puts (  "\nEnter  number  of  sorties  desired  for  Task  3"  ); 
scanf (  " %d"  ,  &Task3  ); 
clear_kb ( ) ; 

/*  Calculate  decision  space  cardinality  using  equation  3.2 

in  thesis  */ 

cardDS  =  factorial (  (double) (  Taskl  +  numMRRs  -  1  )  )  /  ( 

factorial (  (double)  (  Taskl  )  )  *  factorial  (  (double)  (  numMRRs  -  1  )  )) 

*  factorial (  (double) (  Task2  +  numMRRs  -  1  )  )  / 

(  factorial (  (double) (  Task2  )  )  *  factorial (  (double) (  numMRRs  -  1  )  )) 

*  factorial (  (double) (  Task3  +  numMRRs  -  1  )  )  / 

(  factorial (  (double) (  Task3  )  )  *  factorial (  (double) (  numMRRs  -  1  ) 

>>  ; 


printf (  "\nCalculated  decision  space  cardinality  =  %f\n\n", 

cardDS  ) ; 

printf (  "Press  'r1  to  retry,  or  press  Enter  to  continue . \n" 

)  ; 

retry  =  getchO; 
if  (  retry  !=  '  r'  )  { 

break; 

}  //  end  if 
}  //  end  while 

clear_kb ( ) ; 

puts ("Enter  filename."); 
gets (  filename  ) ; 

puts ( " \nEnter  progress  report  increment,  e.g.  100  for  every  100th 
solution . \n" ) ; 

scanf  (  "%d" ,  &t  ) ; 

puts ( "\nEnter  's'  for  screen  output  of  progress,  or  press  Enter  to 
continue . \n" ) ; 

output  =  getchO; 

/*  Open  progress  file  */ 

if  (  (  fp2  =  fopen(  progressf ile ,  "w"  )  )  ==  NULL  )  { 

fprintf (  stderr,  "Error  opening  progress  file."  ); 
exit ( 1) ; 

}  //  end  if 
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/*  Write  header  to  progress  file  */ 

fprintf (  fp2,  "Solution  number\tElapsed  time  (seconds) \n"  ); 
f close  (  fp2  )  ; 

/*  allocate  memory  for  1st  nondominated  solution  (3  fitness  values 
4-  15  DVs )  */ 

Pcurrent  =  (double  *)  realloc (NULL,  (18  *  sizeof (double) )) ; 

/*  memory  allocation  test  */ 
if  (  Pcurrent  ==  NULL  )  { 

puts ( "Memory  allocation  error.") ; 
exit (1) ; 

}  //  end  if 

/*  initialize  1st  solution  with  negatives  if  maximizing  */ 

Pcurrent [0]  =  -1;  //  max 
for  (  i  =  1;  i  <=  17;  i++  )  { 

Pcurrent  [i]  =  100000000;  //  min 


}  //  end 

for 

sizePcurrent  = 

maxXll 

— 

Taskl ; 

maxX12 

= 

Taskl ; 

maxX13 

= 

Taskl ; 

maxX14 

= 

Taskl ; 

maxXIS 

= 

Taskl ; 

maxX21 

= 

Task2 ; 

maxX22 

= 

Task2 ; 

maxX23 

= 

Task2 ; 

maxX24 

= 

Task2 ; 

maxX25 

= 

Task2 ; 

maxX3 1 

= 

Task3 ; 

maxX32 

= 

Task3 ; 

maxX33 

= 

Task3 ; 

maxX34 

= 

Task3 ; 

maxX35 

= 

Task3 ; 

/*  Record  start  time  */ 
start  =  time ( 0 ) ; 

/*  increment  1st  DV  */ 

for  (il  =  0;  il  <=  maxXll;  il++)  { 

/*  increment  2nd  DV  */ 

for  ( i 2  =  0 ;  i2  <=  maxX12 ;  i2++)  { 

/*  increment  3rd  DV  */ 

for  (i3  =  0;  i3  <=  maxX13;  i3++)  { 

/*  increment  4th  DV  */ 

for  ( i4  =0;  i4  <=  maxX14 ;  i4++)  { 

/*  increment  5th  DV  */ 

for  ( i 5  =  0;  i5  <=  maxX15;  i5++)  { 
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/*  increment  6th  DV  */ 

for  (i6  =  0;  i6  <=  maxX21;  i6++)  { 

/*  increment  7th  DV  */ 

for  (i7  =0;  i7  <=  maxX22;  i7++)  { 

/*  increment  8th  DV  */ 

for  (i8  =  0;  i8  <  =  maxX23;  i8++)  { 

/*  increment  9th  DV  */ 

for  (i9  =  0 ;  i9  <=  maxX24;  i9++)  { 

/*  increment  10th  DV  */ 

for  (ilO  =  0;  ilO  <=  maxX25;  il0++)  { 

/*  increment  11th  DV  */ 

for  (ill  =  0;  ill  <=  maxX31;  ill++)  { 

/*  increment  12th  DV  */ 

for  (il2  =  0/  il2  <=  maxX32;  il2++)  { 

/*  increment  13th  DV  */ 

for  (il3  =  0/  il3  <=  maxX33;  il3++)  { 

/*  increment  14th  DV  */ 

for  (il4  =  0;  il4  <=  maxX34;  il4++)  { 

/*  increment  15th  DV  */ 

for  (il5  =  0;  il5  <=  maxX35;  il5++)  { 

/*  test  that  sum  of  like-task  DVs  is  exactly  the  DV's 

task  preference  */ 

if  (  il  +  i2  +  i3  +  i4  +  i5  ==  Taskl 

ScSc  i6  +  i7  +  i8  +  i9  +  ilO  ==  Task2 
&&  ill  +  il2  +  il3  +il4  +  il5  ==  Task3)  { 

/*  evaluate  objective  functions  */ 
evaluate (  fitness,  il,  i2,  i3,  i4,  i5,  i6,  i7, 
i8 ,  i9 ,  ilO,  ill,  il2,  il3,  il4,  il5  ); 

/*  initialize  counter  */ 
n  =  0  ; 

/*  initialize  nondomination  flag;  it's  only 
changed  if  new  solution  is  dominated  */ 

nondominated  =  1; 

/*  while  n  <  the  number  of  solutions  in  Pcurrent 

*/ 

while  (  n  <  sizePcurrent  )  { 

z  =  pareto (  n  ) ; 
switch  (  z  )  { 
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/*  Solution  n  in  Pcurrent  is  dominated  by 
new  solution  &  flagged  for  removal  by  setting  the  1st 

fitness  value  to  a  negative  number 

*/ 

case  1 :  { 

Pcurrent [18  *  n]  =  -1; 
break; 

}  //  end  case  1 

/*  new  solution  is  dominated,  stop  Pareto 
testing,  set  domination  flag  */ 

case  2 :  { 

n  =  (sizePcurrent ) ; 
nondominated  =  0; 
break; 

}  //  end  case  2 

/*  new  solution  is  indifferent  and  will  be 

added  to  Pcurrent  */ 

default:  { 

break; 

}  //  end  default 
}  //  end  switch 
n++  ; 

}  //  end  while 

/*  increment  count  of  solutions  evaluated  */ 
countDS ++ ; 


// 


*/ 


solutions  */ 


==  NULL  )  { 

progress  file .  " 


/*  Diagnostic  line  */ 

printf (  "%f  solutions  evaluated\n" ,  countDS  ); 

/*  print  progress  every  t  solutions  evaluated  */ 
if  (  mod(  countDS, t  )  ==  0  )  { 

/*  record  time  every  t  solutions  */ 
time2t  =  time(0); 

/*  Calculate  time  to  process  100  solutions 

duration  =  difftime(  time2t,  start  ); 

/*  Calculate  percent  complete  of  all 

progress  =  100  *  (  countDS  /  cardDS  ) ; 

/*  Reopen  progress  file  for  appending  */ 
if  (  (  fp2  =  fopen(  progressf ile ,  "a"  )  ) 

fprintf (  stderr,  "Error  opening 

exit  (1) ; 

}  //  end  if 

/*  Write  progress  to  progress  file  */ 
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duration  ) ; 


desired  */ 

seconds . \nH , 


fprintf (  fp2,  "%f\t%f\n";  countDS, 
f close (  fp2  ) ; 

/*  Used  when  screen  output  of  progress  is 

if  (  output  ==  's'  )  { 

printf ( "Time  to  solution  %f  =  %f 

countDS,  duration  ) ; 

printf (  "%f  percent  completed\n\n" , 


progress 


}  //  end  if 
}  //  end  if 


/*  if  new  solution  is  nondominated,  add  new 
solution  to  the  end  of  Pcurrent  */ 

if  (  nondominated  ==  1  )  { 


/*  here's  some  more  memory;  increases  by  3 

doubles  worth  of  bytes  */ 

Pcurrent  =  (double  *)  realloc (Pcurrent , 

(18  *  sizePcurrent  *  sizeof (double)  +  18  *  sizeof (double) )) ; 


allocated  */ 


fitness [1] ; 
fitness [2] ; 


Pcurrent  */ 


/*  memory  allocation  test  */ 
if  (  Pcurrent  ==  NULL  )  { 

puts ("Memory  allocation  error.") ; 
exit (1) ; 

}  //  end  if 

/*  set  values  for  new  memory  that  was 
Pcurrent [18  *  sizePcurrent]  =  f itness  [0] ; 


Pcurrent [18 

* 

sizePcurrent 

+ 

i] 

= 

Pcurrent [18 

* 

sizePcurrent 

+ 

2] 

= 

Pcurrent [18 

* 

sizePcurrent 

+ 

3] 

=  il 

Pcurrent  [18 

* 

sizePcurrent 

+ 

4] 

=  i2 

Pcurrent [18 

* 

sizePcurrent 

+ 

5] 

=  i3, 

Pcurrent [18 

* 

sizePcurrent 

+ 

6] 

=  i4 

Pcurrent [18 

* 

sizePcurrent 

+ 

7] 

= 

Pcurrent [18 

* 

sizePcurrent 

+ 

8] 

=  i6 

Pcurrent [18 

★ 

sizePcurrent 

+ 

9] 

=  i7 ; 

Pcurrent [18 

* 

sizePcurrent 

+ 

10] 

=  i8  ; 

Pcurrent [18 

* 

sizePcurrent 

+ 

11] 

=  i9; 

Pcurrent [18 

★ 

sizePcurrent 

+ 

12] 

=  ilO 

Pcurrent [18 

* 

sizePcurrent 

+ 

13] 

=  ill 

Pcurrent [18 

* 

sizePcurrent 

+ 

14] 

=  il2 

Pcurrent [18 

★ 

sizePcurrent 

+ 

15] 

=  il3 

Pcurrent [18 

* 

sizePcurrent 

+ 

16] 

=  il4 

Pcurrent [18 

* 

sizePcurrent 

+ 

17] 

=  i!5 

/*  increment  count  of  solutions  in 


sizePcurrent++ ; 
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}  //  end  if 
}  //  end  if 


} 

// 

end 

for 

} 

// 

end 

for 

} 

// 

end 

for 

} 

// 

end 

for 

} 

// 

end 

for 

} 

// 

end 

for 

} 

// 

end 

for 

} 

// 

end 

for 

} 

// 

end 

for 

} 

// 

end 

for 

} 

// 

end 

for 

} 

// 

end 

for 

} 

// 

end 

for 

} 

// 

end 

for 

}  //  end  for 

/*  record  finish  time  */ 
finish  =  time (0)  ; 

/*  calculate  total  execution  time  */ 
duration  =  diff time ( finish,  start); 

/*  Reopen  progress  file  for  appending  */ 

if  (  {  fp2  =  fopen(  progressf ile ,  "a"  )  )  ==  NULL  )  { 

fprintf (  stderr,  "Error  opening  progress  file."  ); 
exit (1)  ; 

}  //  end  if 

/*  Write  to  progress  file  */ 

fprintf  (  fp2,  "%f\t%f\n,,/  countDS,  duration  ); 
f close  (  fp2  )  ; 

/*  Used  if  screen  output  of  progress  is  desired  */ 
if  (  output  ==  ,sI  )  { 

printf (  "\nProgram  execution  time  =  %f  seconds. \n",  duration 

)  ; 

printf (  "Calculated  decision  space  cardinality  =  %f.\n", 

cardDS  ) ; 

printf (  "Decision  space  cardinality  =  %f.\n"/  countDS  ); 

}  //  end  if 

/*  Open  report  file  */ 

if  (  (  fpl  =  fopen(  filename,  "w"  )  )  ==  NULL  )  { 

fprintf (stderr,  "Error  opening  file  %s.",  filename); 
exit (1) ; 

}  //  end  if 

/*  Write  to  report  file  */ 

fprintf (  fpl, "OBJECTIVE  l\tOBJECTIVE  2\tOBJECTIVE 
3\txll\txl2\txl3\txl4\txl5\tx2l\tx22\tx23\tx24\tx25\tx3l\tx32\tx33\tx34\ 
tx35\n") ; 

for  (  i  =  0;  i  <  sizePcurrent ;  i++  )  { 

if  (  Pcurrent [18  *  i]  >=  0  ) 
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fprintf ( 

fpl, "%lf \t%lf \t%f\t%f \t%f \t%f\t%f \t%f\t%f \t%f \t%f\t%f \t%f \t%f \t%f \t%f \t% 
f \t%f \n" ,  Pcurrent [18  *  i]  ,  Pcurrent [18  *  i  +  1] ,  Pcurrent [18  *  i  +  2]  , 
Pcurrent[18  *  i  +  3] ,  Pcurrent [18  *  i  +  4] ,  Pcurrent  [18  *  i  +  5] , 

Pcurrent  [18  *  i  +  6],  Pcurrent [18  *  i  +  7] ,  Pcurrent [18  *  i  +  8] , 

Pcurrent [18  *  i  +  9] ,  Pcurrent [18  *  i  +  10],  Pcurrent  [18  *  i  +  11], 
Pcurrent [18  *  i  +  12],  Pcurrent [18  *  i  +  13],  Pcurrent [18  *  i  +  14], 

Pcurrent [18  *  i  +  15],  Pcurrent [18  *  i  +  16],  Pcurrent [18  *  i  +  17]  ); 

}  //  end  for 

/*  close  files  and  release  allocated  memory  */ 
f close (  fpl  ) ; 
free (  Pcurrent  )  ; 

}  //  end  main 


void  evaluate (double  f [] ,  int  il,  int  i2 ,  int  i3,  int  i4,  int  i5,  int 
i6,  int  i7 ,  int  i8,  int  i9,  int  ilO,  int  ill,  int  il2,  int  il3,  int  il4, 
int  il5  ) 

{ 

f  [0]  =  . 8*il  +  .  3*i2  +  . 6*i3  +  ,001*i4  +  ,001*i5  +  .4*i6  +  .8*i7  + 

. 6*i8  +  . 001*i 9  +  . 001*il0  +  .001*ill  +  ,001*il2  +  .l*il3  +  .8*il4  + 

.  4  *  i  1 5 ; 

f[l]  =  20.2*  (il  +  i6  +  ill)  +  28.5*(i2  +  i7  +  il2)  +  35.7*(i3  +  i8 

+  il3)  +  19.9*  (i4  +  i9  +  il4)  +  22.5*(i5  +  ilO  +  il5); 

f  [2]  =  1650*  (il  +  i6  +  ill)  +  2475*  ( i 2  +  i7  +  il2)  +  2887.5*(i3  + 

i8  +  il3 )  +  1705*  (i4  +  i9  +  il4)  +  2200* (i5  +  ilO  +  il5) ; 

} 


/*  Pareto  dominance  truth  table  for  maximizing  objective  1,  minimize 
objectives  2  and  3  */ 
int  pareto (int  n) 

{ 

int  m; 


if  ( 
fitness  [1] 

li 

fitness  [1] 
fitness  [1] 
fitness  [1] 
fitness  [1] 
fitness  [1] 
fitness  [1] 

}  // 


Pcurrent [18*n] 
&&  Pcurrent  [18 
Pcurrent [18*n] 
ScSc  Pcurrent  [18 
Pcurrent [18*n] 
&&  Pcurrent [18 
Pcurrent [18*n] 
ScSc  Pcurrent  [18 
Pcurrent [18*n] 
&&  Pcurrent [18 
Pcurrent [18*n] 
ScSc  Pcurrent  [18 
Pcurrent [18*n] 
ScSc  Pcurrent  [18 
/*  solution  n 


==  fitness  [0]  ScSc  Pcurrent  [18*n  + 
*n  +  2]  >  fitness [2] 

==  fitness  [0]  ScSc  Pcurrent  [18*n  + 
*n  +  2]  ==  fitness  [2] 

==  fitness  [0]  ScSc  Pcurrent  [18*n  + 
*n  +  2]  >  fitness  [2] 

<  fitness  [0]  ScSc  Pcurrent 
*n  +  2]  ==  fitness [2] 

<  fitness  [0]  ScSc  Pcurrent 

*n  +  2]  >  fitness  [2] 

<  fitness  [0]  ScSc  Pcurrent 
*n  +  2]  ==  fitness  [2] 

<  fitness  [0]  ScSc  Pcurrent 

*n  +  2]  >  fitness [2]  )  { 

in  Pcurrent  is  dominated 


1] 


1] 


1]  == 


[18*n 

+ 

i] 

== 

[18*n 

+ 

i] 

> 

[18*n 

+ 

i] 

> 

[18*n 

+ 

i] 

=  = 

by  new 

solution  */ 

m  =  1  ; 
return  m; 
end  if 
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if  ( 
fitness  [1] 

fitness  [1] 

II 

fitness  [1] 

ii 

fitness  [1] 
fitness  [1] 
fitness  [1] 
fitness  [1] 

}  // 


Pcurrent [18*n]  >  fitness [0]  &&  Pcurrent  [18*n  +  1]  == 

ScSc  Pcurrent  [18*n  +  2]  ==  fitness  [2] 

Pcurrent [18*n]  ==  fitness [0]  ScSc  Pcurrent [18*n  +  1]  == 

&&  Pcurrent  [18*n  +  2]  <  fitness [2] 

Pcurrent [18*n]  ==  fitness [0]  &&  Pcurrent [18*n  +  1]  < 

ScSc  Pcurrent  [18*n  +  2]  ==  fitness  [2] 

Pcurrent [18*n]  ==  fitness [0]  &&  Pcurrent [18*n  +  1]  < 

ScSc  Pcurrent  [18*n  +  2]  <  fitness  [2] 

Pcurrent  [18*n]  >  fitness  [0]  ScSc  Pcurrent  [18 *n  +  1]  == 

ScSc  Pcurrent  [18*n  +  2]  <  fitness  [2] 

Pcurrent [18*n]  >  fitness [0]  &&  Pcurrent  [18*n  +  1]  < 

ScSc  Pcurrent  [18*n  +  2]  ==  fitness  [2] 

Pcurrent  [18*n]  >  fitness  [0]  ScSc  Pcurrent  [18*n  +  1]  < 

ScSc  Pcurrent  [18*n  +  2]  <  fitness  [2]  )  { 

/*  new  solution  is  dominated  by  solution  n  in  Pcurrent 
m  =  2  ; 
return  m; 
end  if 


*/ 


} 


m  =  3  ; 
return  re¬ 


double  factorial (  double  a  ) 

{ 

if  (  a  ==  1  ) 

return  1; 

else  { 

a  *=  factorial (  a  -  1  ) ; 
return  a; 

} 

} 

double  mod(  double  a,  int  b  ) 

{ 

double  c; 

/*  using  ceil  function  since  I  can't  find  a  rounding  funtion  */ 
c  =  ceil  (  f mod (  a,  (double) b  )  *  b  ); 

return  c; 

} 

void  clear__kb(  void  ) 

{ 

/*  Clears  stdin  of  any  waiting  characters.  */ 
char  junk [80] ; 
gets (  junk  ); 

} 
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Appendix  D:  Ptrue  and  PFtrue  for  Table  4,  Index  1 
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Appendix  E:  Source  Code  for  Pareto_processing.c 


Pareto_processing.c  V2.2 
2/18/01 

Dave  Wakefield  &  Jesse  Zydallis.  Thanks  to  Matt  Johnson. 

Side  constraint,  Pareto  check,  &  clone  check  program. 

This  program  will  take  the  data  points  and  remove  those  that  are 
1)  infeasible,  2)  dominated,  and  3)  clones.  These  three  routines  are 
run  successively  on  an  input  file,  generating  a  report  after  each 
routine.  However,  input  and  output  file  names  are  hard  coded,  along 
with  defined  parameters.  The  data  array  size  in  the  main  function  is 
also  hard  coded. 

=  =  =  =  „«====  =======«===  =  */ 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  <string.h> 

#def ine  NUM_DVS  15; 

#def ine  NUM_FUNCS  3; 

#def ine  MAX_COLS  18; 

#define  MAX_PTS  40000;  /*  set  very  high  since  number  of  solutions  is 
unknown  * / 

#def ine  PAGEWIDTH  63; 

#def ine  REPS  30; 

#def ine  MAXIMIZATION_FLAG  0;  /*  1  means  max  */ 

int  num_DVs,  num_funcs; 
num___DVs  =  NUM_DVS  ; 
num__f  uncs  =  NUM_FUNCS ; 

int  constraint (  double  *ex  ) ; 

void  is_par(  char  ^filename,  int  rows,  int  file_num  ); 
void  is__clone  (  char  ^filename,  int  rows,  int  file_num  ); 

main()  { 

register  int  i,  j,  m;  /*  for  loop  indices  */ 

int  pass;  /*  represents  boolean  result  from  constraint  test  */ 
int  nind,  column;  /*  counts  input  file  rows  and  columns  */ 
int  total;  /*  used  to  detect  end  of  input  file  row  condition  */ 
int  z;  /*  counts  number  of  feasible  solutions  */ 
int  reps ,  max_pts ; 

double  *ex;  /*  array  holding  DVs  for  single  solution  */ 
double  number;  /*  temp  var  used  to  read  in  data  from  input  file 

*/ 

FILE  *fp,  *fp2 ;  /*  fp  reads,  fp2  writes  */ 

char  filel[30],  file2  [30]  ;  /*  filel  and  file2  names  */ 

double  data[40000]  [18];  /*  array  to  hold  input  file  data  */ 

total  =  MAX__C0LS; 
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reps  =  REPS; 
max_jpts  =  MAX_PTS; 


/*  outer  loop  used  to  process  input  files  from  all  replications 

*/ 

for  (  m  =  1;  m  <=  reps;  m++  )  { 

/*  set  filel  to  desired  input  file  name  and  open  it  for 

reading  */ 

if  (  m  <  10  )  { 

sprintf  (  filel,  "el_r0%d . pts " ,  m  ); 

}  else  { 

sprintf  (  filel,  "el__r%d . pts "  ,  m  ); 

}  /*  end  if  */ 

if  (  (  fp  =  fopen(  filel,  "r"))  =  =  NULL  )  { 

fprintf (  stderr,  "Error  opening  file  %s . " ,  filel  ); 
exit (1) ; 

}  /*  end  if  */ 

column  =  nind  =  z  =  0; 


*/ 


/*  read  in  data  from  input  file  */ 
while ( (fscanf (fp,  "%lf " ,  ^number) ) !=EOF) { 
data [nind] [column]  =  number; 
column++ ; 

column%=total;  /*  detects  end  of  row  */ 
if  (column==0 )  {  /*  if  end  of  row,  move  to  next  row 

nind++ ; 


if  (nind==max_pts)  { 

print f ( "Too  many  data  points  for 

side_constraint . \n" ) ; 

printf ( "Increase  1 MAX_PTS ' \nn ) ; 
exit (1) ; 

}  /*  end  if  */ 

}  /*  end  if  */ 

}  /*  end  while  */ 
f close  (fp)  ; 

/*  allocate  memory  to  hold  a  row  of  DVs  and  OFs  */ 
if ( ! (  ex  =  (double  *)malloc(  total*sizeof (double)  )  ))  { 

fprintf (  stderr,  "Insufficient  memory  for  variable, 

ex . \n" ) ; 

}  /*  end  if  */ 

/*  set  f ile2  to  desired  out  file  name  and  open  it  for 

writing  */ 

if  (  m  <  10  )  { 

sprintf  (  file2,  "el__rO%d .  f eas  "  ,  m  ); 

}  else  { 

sprintf (  file2,  "el_r%d . f eas" ,  m  ); 

}  /*  end  if  */ 

if  ((  fp2  =  fopen(  file2,  "w" ))  ==  NULL  )  { 

fprintf (  stderr,  "Error  opening  file  %s.",  file2  ); 
exit (1) ; 
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}  /*  end  if  */ 

/*  this  loop  checks  each  row  of  data  array  for  feasibility 

*/ 

for  (  i  =  0;  i  <  nind;  i++  )  { 

for  (  j  =  0;  j  <  num_DVs;  j++  )  { 

*  (ex+j  )  =  data  [i]  [j  ]  ; 

}  /*  end  for  */ 

pass  =  constraint (ex) ; 
if  (  pass  ==  1  )  { 

for  (  j  =  0 ;  j<  total;  j++  )  { 

fprintf (  fp2,  »%3.14f  %  data[i][j]  ); 

}  /*  end  for  */ 
fprintf (  fp2,  M\nn  ); 
z+  +  ; 

}  /*  end  if  */ 

}  /*  end  for  */ 
f close  (  fp2  ) ; 

/*  pass  feasible  points  file  to  Pareto  check  */ 
is_par(  file2,  z,  m  ); 

}  /*  end  for  */ 

}  /*  end  main  */ 

int  constraint (double  *ex) 

{ 

/*  print f ( "ex  =  %f\n",*ex);  */ 

/♦insert  function  here*/ 

if  (  ( (*ex  +  *(ex+l)  +  *(ex+2)  +  *(ex+3)  +  *(ex+4))  ==  10.0)  && 

((*(ex+5)  +  *(ex+6)  +  * (ex+7 )  +  * (ex+8)  +  * (ex+9) )  ==  5.0) 

ScSc 

( ( * (ex+10 )  +  *(ex+ll)  +  *(ex+12)  +  *(ex+13)  +  *  (ex+14) )  == 

1.0) 

)  { 

return  (1) ; 

}  else  { 

return  (0) ; 

/*  if  (  (  *ex* (*ex) +  * (ex+1) * (* (ex+1) ) <=225)  &&  (  (  *ex-3* (* (ex+1) ) 

) <= - 10  )  )*/ 

} 

function  :  is_par 

purpose  :  finds  pareto  optimal  points 
developed  :  2001  from  newis_par.c  by  Matt  Johnson 
modified  by  Dave  Wakefield  &  Jesse  Zydallis 
void  is_par (  char  ^filename,  int  rows,  int  file_num  ) 
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register  int  i,  j,  k; 

int  flagl,  flag2,  count,  rank,  total,  maximizationf lag ; 

FILE  *fp,  *fp2 ; 

double  number,  *answer_data ,  *answer_ptr; 
double  *data,  *ptrl,  *ptr2,  *ptr3; 
double  frac_par; 
char  filel [30] ; 

maximizationflag  «  MAX I M I Z AT I ON_ F LAG ; 
total  =  num_DVs  +  num^funcs; 

printf ("Inside  Pareto  Analysis  Routine !!! \n" ) ; 

data  =  (double  *) malloc (rows*total*sizeof (double)  )  ; 

if  (data==NULL)  { 

fprintf (  stderr,"Not  enough  memory  for  output  data.\n"); 

}  /*  end  if  */ 

answer^data  =  (double  *) malloc (rows*total*sizeof (double)  )  ; 

if  (answer_data==NULL)  { 

fprintf (  stderr, "Not  enough  memory  for  output  data.\n"); 

}  /*  end  if  */ 

/*  open  feasible  data  file,  *.feas  */ 

if  ( (  fp  =  fopen(  filename,  "r"))  ==  NULL  )  { 

fprintf (  stderr,  "Error  opening  file  %s . " ,  filename  ); 
exit (1) ; 

}  /*  end  if  */ 

ptrl  ==  data; 

while ( (fscanf (fp,  "%lf",  &number) ) ! =EOF) { 

*ptrl=number ;  /*input  all  the  data*/ 

ptrl++; 

}  /*  end  while  */ 
f close (fp) ; 

count =0 ;  /*will  be  the  number  of  pareto  optimal 

pts  */ 

if  (maximizationf lag==0)  {  /*a  minimization  problem*/ 

for  ( i  =  0 ;  i<rows;  i++)  { 

ptrl=data+i*total;  /*beginning  of  each  row*/ 

rank=rows ; 

ptr2  =ptrl+num_DVs ;  /*start  at  fl  for  each  row*/ 

for  ( j  =  0 ;  j  <rows ;  j++)  { 

ptr3=data+j*total+num_DVs;  /*go  through  all 

rows*/ 

f lagl=f lag2=0 ; 

for  (k=0 ; k  <  num_funcs;  k++)  {  /*go 

through  all  the  functions  of  a  row*/ 

if  ( (* (ptr2+k) ) < (* (ptr3+k) ) )  { 
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flagl++;  /*row  can't  be 

dominated  so  break*/ 

break; 

}  /*  end  if  */ 

if  ( (* (ptr2+k) ) == (* (ptr3+k) ) )  { 
f lag2++; 

}  /*  end  if  */ 

}  /*  end  for  */ 

if  (flagl>0  | |  flag2  ==  num_funcs)  {  /*non- 

dominated  | |  same  row*/ 

rank--; 

}  else  { 

break;  /*a  row  is  dominated, 

move  on*/ 

}  /*  end  if  */ 

}  /*  end  for  */ 
if  (rank==0)  { 

answer__ptr=answer_data+count  *  total ;  /*move  down 

count  rows*/ 

for(j=0;  jctotal;  j++)  { 

* (answer_ptr+j ) =* (ptrl-f-j ) ;  /*take 

points  and  function  values*/ 

}  /*  end  for  */ 
count++ ; 

}  /*  end  if  */ 

}  /*  end  for  */ 

}  else  {  /*a  maximization  problem*/ 
for  (i=0;  i<rows;  i++) 

{ 

ptrl=data+i*total ;  / *beginning  of  each  row*/ 

rank=rows ; 

ptr2  =  ptrl  +  num_DVs;  /*  start  at  fl  for  each  row 


7 


rows*/ 


for  ( j  =  0 ;  j<rows;  j+  +  ) 

{ 

ptr3=data+j *total  +  num_DVs;  /*go  through  all 


f lagl=f lag2=0 ; 

for  (k=0 ; k  <  num_funcs;  k++)  /*go  through  all 


the  functions  of  a  row*/ 


{ 


dominated  so  break*/ 


dominated  | |  same  row*/ 


if  ( (* (ptr2+k) ) > (* (ptr3+k) ) ) 

{ 

flagl++;  /*row  can't  be 
break; 

} 

if  ( (* (ptr2+k) ) == (* (ptr3+k) ) ) 

{ 

f lag2++ ; 

} 

} 

if  (flagl>0  ||  flag2  ==  num__funcs  )/*non- 


{ 
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rank- - ; 


} 

else 

{ 

break;  /*a  row  is  dominated, 


if  (rank==0) 

{ 

answer__ptr=answer__data+count  *  total ;  /*move  down 

count  rows*/ 

for(j=0;  jctotal;  j++) 

{ 

* (answer_ptr+j ) =* (ptrl+j ) ;  /*take 

points  and  function  values*/ 

} 

count++; 

} 

} 

} 

free (data) ; 

if  (  filejnum  <  10  )  { 

sprintf(  filel,  "el_r0%d.prto" ,  file_num  ); 

}  else  { 

sprintf (  filel,  "el_r%d .prto" ,  file_num  ); 

}  /*  end  if  */ 

if  ((  fp2  =  f open (  filel,  "w"))  ==  NULL  )  { 

fprintf (  stderr,  "Error  opening  file  %s." ,  filel  ); 
exit (1) ; 

}  /*  end  if  */ 

for  ( i = 0 ;  i< (total*count) ;  i++)  {  /*  write  the  data  to  the 

*.prto  file  */ 

fprintf (fp2,  "%3.14f  ",  * (answer_data+i ) ) ; 

if  (  (i+1) %total==0  )  {  /*would  mean  finished  a  whole  row*/ 

fprintf (fp2 , "\n") ; 

}  /*  end  if  */ 

}  /*  end  for  */ 
f close ( fp2 ) ; 

free (answer_data) ;  /*free  memory*/ 

f  rac_jpar=  (double)  count/  (double)  rows  ; 

printf("%s  fraction  of  pareto  optimal  points  =  %3.14f\n" ,  filel, 
frac_par) ; 

/*  pass  Pareto  points  to  is_clone  to  remove  duplicate  points  */ 
is_clone(  filel,  count,  file_num  ); 

}  /*  end  is_par  */ 

function  :  is_clone 

purpose  :  finds  duplicate  points 
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developed  :  2001  from  newis_par.c  by  Matt  Johnson 
modified  by  Dave  Wakefield  &  Jesse  Zydallis 
void  is_clone(  char  ^filename,  int  rows,  int  file_num  )  { 

register  int  i,  j,  k; 

int  total;  /*  number  of  elements  in  a  row  */ 

int  count;  /*  counts  the  number  of  unique  pareto  optimal  pts  */ 
FILE  *fp,  *fp2; 

double  number,  *answer_data ,  *answer__ptr ; 
double  *data,  *ptrl,  *ptr2; 
char  f ilel  [30]  ; 

total  =  num_DVs  +  num_funcs; 

printf (  "Inside  No  Clone  Routine! ! !\n"  ); 

data  =  (double  *)malloc(  rows  *  total  *  sizeof (double)  )  ; 

if  (  data  ==  NULL  )  { 

fprintf (  stderr, "Not  enough  memory  for  output  data.\n"  ); 

}  /*  end  if  */ 

answer_data  =  (double  *)malloc(  rows  *  total  *  sizeof (double)  )  ; 

if  (  answer__data  ==  NULL  )  { 

fprintf (  stderr, "Not  enough  memory  for  output  data.\n"  ); 

}  /*  end  if  */ 

/*  open  Pareto  points  data  file,  *.prto  */ 
if  (  (  fp  =  fopen(  filename,  "r"))  ==  NULL  )  { 

fprintf (  stderr,  "Error  opening  file  %s.",  filename  ); 
exit (1) ; 

}  /*  end  if  */ 

ptrl  =  data; 

while (  (  f scanf ( fp ,  "%lf",  ^number)  )  !=  EOF  ){ 

*ptrl  =  number;  /*  input  all  the  data  */ 
ptrl++; 

}  /*  end  while  */ 
f close (  fp  ) ; 
count  =  0 ; 

/*  next  3  for  loops  compare  COLUMN  k  of  ROW  i  against  COLUMN  k  of 
ROW  j  */ 

for  (  i  =  0;  i  <  rows;  i++  )  { 

ptrl  =  data  +  i  *  total;  /*  point  to  beginning  of  ROW  i  */ 
if  (  *ptrl  !=-l)  {  /*  -1  means  the  row  is  a  clone  and  it 

won't  be  tested  */ 

for  (  j  =0;  j  <  rows;  j++  )  { 

ptr2  =  data  +  j  *  total;  /*  point  to  beginning 

of  ROW  j  */ 
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if  (  i  !=  j 

compare  a  row  against  itself;  don't 

for  ( 

k)  )  )  {  /*  point  to  COLUMN  k  */ 

j ,  so  break  to  next  j  */ 


&&  *ptr2  ! =  -1  )  {  /*  don't  want  to 

want  a  duplicate  row  */ 
k  =  0;  k  <  num_DVs;  k++  )  { 

if  (  (  * (ptrl  +  k)  )  !=  (  * (ptr2  + 

/*  ROW  i  is  not  a  duplicate  of  ROW 

break; 

}  /*  end  if  */ 


if  (  k  ==  (  num_DVs  - 

last  DV  has  been  checked,  ROW  j  is  a  clone  of  ROW  i  */ 

data[j  *  total] 


ROW  j 


clone  flag  */ 

} 


}  /*  end  if  */ 
}  /*  end  for  */ 

}  /*  end  if  */ 

/*  end  for  */ 


1  ))  {  /*  if 

=  -1;  /*  set 


/*  if  you  get  this  far,  then  ROW  i,  having  flagged 
any  clones  of  itself,  is  unique,  so  copy  it  to  answer_data  */ 

answer_ptr  =  answer_data  +  count  *  total;  /*  move 
down  <count>  rows  to  append  ROW  i  to  answer__data  */ 

for  (  j  =  0;  j  <  total;  j++  )  { 

* (answer_ptr+j )  =  * (ptrl+j ) ;  /*  copy  DVs 

and  function  values  */ 

}  /*  end  for  */ 
count++; 

}  /*  end  if  */ 

}  /*  end  for  */ 


free (  data  ) ; 


if  (  file_num  <  10  )  { 

sprintf (  filel,  "el_r0%d . front " ,  file_num  ); 

}  else  { 

sprintf (  filel,  "el_r%d. front" ,  file_num  ); 

}  /*  end  if  */ 

if  ((  fp2  =  fopen(  filel,  "w") )  ==  NULL  )  { 

fprintf (  stderr,  "Error  opening  file  %s.",  filel  ); 
exit ( 1) ; 

}  /*  end  if  */ 

for  (i=0;  i<  (total*count)  ;  i++)  {  /*  write  data  to  the 

*. front  file  */ 

fprintf (fp2,  "%3.14f  ",  *  (answer_data+i) ) ; 

if  (  (i+1) %total==0  )  {  /*would  mean  finished  a  whole  row*/ 

fprintf (fp2 , "\n" ) ; 

}  /*  end  if  */ 

}  /*  end  for  */ 
f close ( fp2 ) ; 

f ree (answer_data) ;  /*free  memory*/ 

printf  (  "%s  non-dominated  set  cardinality  =  %d\n" ,  filel,  count 

)  ; 

}  /*  end  is_clone  */ 
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Appendix  F:  Raw  Data  and  Experimental  Statistics 

Table  9.  Raw  Data  for  Final  Generational  Distance 
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Table  10.  Descriptive  Statistics  for  Final  Generational  Distance 
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Count  30  30!  30:  30 

Largest(l)  1590.664698 . 568.9250723  1138.672605:  895.1914708 

Smallest(l)  55.83325178  55.83325178  55.83325178  45.52391375 

Confidence  Level(95.0%)  142.1772638  54.41397388  105.0448362  9  293618901 


Table  1 1 .  Raw  Data  for  Overall  Nondominated  Vector  Generation 
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Table  12.  Descriptive  Statistics  for  Overall  Nondominated  Vector  Generation 
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Confidence  Level(95.0%)  0.756302151  0.8601 11911  0.653415108  0.625853467 


Oneway  Analysis  of  Generational  Distance  By  Parameter 


Figure  13.  Kruskal-Wallis  H- Test  Results  for  Final  Generational  Distance 
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Appendix  G:  3D  Plots  of  PFtrue  and  PFkn0Wn  Using  Alternative  Parameter  Values 


The  following  figures  plot  in  three  dimensions  the  PF,rue  with  the  PFknown 
generated  by  the  implicitly  constraining  MOMGA-II  using  the  alternative  parameter 
values  specified  in  the  experimental  design  section  of  Chapter  III.  The  plots  are  intended 
to  show  that  without  employing  any  explicit  constraint  handling  methods,  the  MOMGA- 
II  output  approximates  the  structure  of  PFm,e. 


Figure  15.  Plot  of  PFtnie  and  PF ' known  for  BB  size  4 
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Figure  1 6.  Plot  of  PFtruc  and  PF known  for  BB  size  8 


Figure  17.  Plot  of  PFtrue  and  PFknow„  for  BB  size  2 
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Figure  18.  Plot  of  PFtrue  and  PFknown  for  Pcut  =  0 


Figure  1 9.  Plot  of  PFtrue  and  PFknown  for  Pcu,  =  0.2 
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Column  4 


Figure  20.  Plot  of  PFtnie  and  PFk  nown  for  P splice  0.8 
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Figure  21 .  Plot  of  PFtnte  and  PF  known  for  P splice  =  0.6 
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Figure  22.  Plot  of  PFtrue  and  PF  known  for  initial  population  size  =  600 
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Figure  23.  Plot  of  PFtrue  and  PFknown  for  initial  population  size  =  1200 
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Appendix  H:  3D  Plots  oiPFknown  for  Resource  Levels  1  -5 


The  following  figures  plot  in  three  dimensions  for  each  Resource  Level  the 
PFknown  generated  by  the  explicitly  constraining  MOMGA-II  using  the  basic  parameter 
values  specified  in  the  experimental  design  section  of  Chapter  III. 
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Figure  24.  PFknown  Plotted  With  PFtnie  for  Resource  Level  1 
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Figure  25.  PF known  Plotted  for  Resource  Level  2 
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Figure  27.  PFknown  Plotted  for  Resource  Level  4 
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Figure  28.  PFknown  Plotted  for  Resource  Level  5 
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